0.25in

# Reflections in excitable media linked to existence and stability of one-
dimensional spiral waves

Stephanie Dodson Department of Mathematics, University of California, Davis,
Davis, CA<EMAIL_ADDRESS>Timothy J. Lewis Department of Mathematics,
University of California, Davis, Davis, CA<EMAIL_ADDRESS>
###### Abstract

When propagated action potentials in cardiac tissue interact with local
heterogeneities, reflected waves can sometimes be induced. These reflected
waves have been associated with the onset of cardiac arrhythmias, and while
their generation is not well understood, their existence is linked to that of
one-dimensional (1D) spiral waves. Thus, understanding the existence and
stability of 1D spirals plays a crucial role in determining the likelihood of
the unwanted reflected pulses. Mathematically, we probe these issues by
viewing the 1D spiral as a time-periodic antisymmetric source defect. Through
a combination of direct numerical simulation and continuation methods, we
investigate existence and stability of a 1D spiral wave in a qualitative ionic
model to determine how the systems propensity for reflections are influenced
by system parameters. Our results support and extend a previous hypothesis
that the 1D spiral is an unstable periodic orbit that emerges through a global
rearrangement of heteroclinic orbits and we identify key parameters and
physiological processes that promote and deter reflection behavior.

Keywords: one-dimensional spiral waves, source-defects, wave reflections,
cardiac arrhythmia, Morris-Lecar

AMS subject classifications: 35B36, 35K57, 35P30, 92-10

##  1 Introduction

The ability of the heart to effectively pump blood to the body relies on the
coordinated spread of electrochemical waves (propagated action potentials)
through the cardiac tissue. If the spread of the waves are disrupted, life-
threatening arrhythmias may arise. The most severe form of arrhythmias involve
reentry behavior, which is composed of self-perpetuating circulating waves
that underlies arrhythmias such as ventricular tachycardia and ventricular
fibrillation. Reentrant arrhythmias are thought to be induced by propagated
action potentials interacting with functional or structural heterogeneities in
the cardiac tissue that lead to symmetry breaking of the wave via local one-
way propagation block.

One of the mechanisms that is thought to induce the onset of reentrant
arrhythmias is reflection [4, 49, 3, 29]. Usually, when a propagated cardiac
action potential interacts with a local tissue heterogeneity, the action
potential is either blocked and annihilated or it successfully propagates
across the heterogeneity. However, in some cases, action potentials can
successfully cross the heterogeneity and then trigger a reflected action
potential, i.e., an action potential that propagates only in the retrograde
direction. A locally reflected pulse breaks the symmetry of the wave of
activation spreading through the heart and thus can initiate reentry.
Reflection has been directly observed in cardiac Purkinje fiber [50, 4, 28,
2], ventricular tissue [5, 44, 17], and atrial tissue [35] in experimental
preparations with regions of depressed excitability or reduced electrical
coupling. More recently, reflection was observed in monolayers of ventricular
cardiomyocytes with sites of abrupt tissue expansions [7]. Action potential
reflection has also been demonstrated in computational models of cardiac
tissue [11, 33, 7, 29]. Furthermore, reflections have been observed in
experimental preparations of neural tissue at the axon level [27, 8, 9] and
the network level [13, 52], as well as corresponding neural models [23, 54,
24].

Previous work began to uncover the dynamical mechanisms underlying reflections
of pulses. In 1996, Ermentrout & Rinzel conducted numerical experiments on a
cable model of an axon using Morris-Lecar dynamics and simulated the effects
of an abrupt change in the diameter on action potential propagation [21]. When
the diameter of the distal portion of the axon was sufficiently large, the
action potential failed to propagate across the heterogeneity; whereas for
sufficiently small distal axon diameters, the action potential successfully
propagated across the heterogeneity. At intermediate diameters, the action
potential propagated across the heterogeneity and then gave rise to a
reflected pulse. Moreover, as the diameter was tuned to be closer and closer
to a critical value, more and more reflections seemed to occur with successive
counter-propagating action potentials. These simulations suggested that an
unstable periodic orbit underlies the reflection behavior, and that the closer
that the incoming wave brings the system to the periodic orbit, the more
reflections occur. The authors pointed out that the unstable periodic orbit
was a so-called one-dimensional (1D) spiral wave that was first described by
Kopell & Howard in [31]. 1D spiral waves are ‘source defects’ [46] that
consist of a non-excited core that sheds anti-phase counter-propagating
pulses.

An example of a series of reflections is given in Figure 1. The pattern has
been dubbed a ‘1D spiral wave’ because its spatial structure resembles a
radial slice through a two-dimensional rotating spiral wave. Waves emerge from
the core on alternating sides and arrange into regular periodic structures as
they propagate into the far-field. However, we emphasize that, despite the
similarities in the names, the 1D and 2D spirals are different mathematical
structures. Our focus here is specifically on the 1D spiral wave.

Figure 1:  (a) Space-time plot of a 1D spiral wave. (b) Example of spatial
profile in 1D spiral. Green arrows above peaks provide direction of motion in
the far-field.

To probe the link between reflection and the 1D spiral, Cytrynbaum & Lewis
[16] considered a homogeneous Morris-Lecar cable model and examined the
interaction of a propagated action potential with a transient refractory
region, i.e., a functional heterogeneity. By tuning the height of the
refractory bump, they observed similar behavior as Ermentrout and Rinzel in
[21]. However, reflection behavior was only seen when the excitability
parameter $\epsilon$ was below a critical value $\epsilon_{*}$. Above this
value $(\epsilon>\epsilon_{*}$), no evidence for the existence of the 1D
spiral was detected, and no signs of reflections were observed (see Figure 4).
Numerical experiments near $\epsilon_{*}$ lead them to postulate that the
unstable periodic orbit is created through a heteroclinic bifurcation at
$\epsilon_{*}$.

Pulse generators that resemble a one-sided 1D spiral have been analyzed in
[40, 38, 48, 53]. These structures were created by a jump-type heterogeneity
and were also found to originate via a heteroclinic bifurcation [53]. However,
there are distinct differences from the scenario of the 1D spiral presented
here. Specifically, the pulse generators of [53] are stable and produce a
variety of heterogeneity-induced patterns that have not been linked to
reflections.

Understanding how system parameters influence the existence of 1D spirals and
the ability of the system to generate reflections is vitally important for
identifying the physiological mechanisms that contribute to reflection-
mediated arrhythmias. Even though the 1D spiral structure is unstable, the
existence of the periodic orbit allows for new spatiotemporal dynamics, in
that initial conditions near it lead to a transient series of spontaneous
reflections. The specific physiological conditions that promote or deter
reflections remain poorly understood, and the dynamical mechanisms creating
the periodic orbit are unverified. It has been suggested [21] that reflection
phenomena occur more robustly when local dynamics exhibit type I excitability
(excitability that is associated with saddle-type threshold behavior) than
type II excitability (quasi-threshold behavior) [43].

Here, we further probe these questions in the Morris-Lecar model. To do so, we
interpret the system’s propensity for reflections in two ways: (1) how
conditions alter the probability of reflections occurring when the 1D spiral
exists, and (2) how conditions modify the 1D spiral existence window. Previous
studies detected the 1D spiral with direct numerical simulations, but these
methods are inefficient for finding unstable solutions. Therefore, to
investigate the 1D spiral existence and properties, we develop more refined
numerical techniques that solve directly for the 1D spiral as a spatiotemporal
equilibrium solution. Numerical methods used and developed will be outlined in
the main text, and full implementation details are included in the Appendix.

This paper is organized as follows. Section 2 describes the Morris-Lecar model
and introduces the mathematical framework of traveling pulses, periodic waves,
reflections, and appearance of 1D spirals. In Section 3, we investigate
scenarios when the 1D spiral does exist and understand how model parameters
influence the likelihood of reflections by linking the stability of the 1D
spiral with the probability of generating reflections. Conditions that alter
the existence of 1D spirals and reflections are probed in Section 4. We
further verify the proposed heteroclinic structure of the 1D spiral
bifurcation by confirming there are no additional local bifurcations of the
slow pulse and accurately computing the period of the orbit near the
bifurcation point. To compute the 1D spiral existence interval and it’s
dependence on parameters, we define novel computational techniques to directly
compute the heteroclinic bifurcation point. Our results are interpreted in the
context of the local excitability type in Section 5 to test the hypothesis
that type I local dynamics promote reflections. Finally, in Section 6 the
stability and existence results are combined with the physiological mechanisms
to discuss which processes could be most involved in generating reflections
and provide predictions on how to decrease reflections in more complex
systems.

##  2 Mathematical Formulation

Our primary goal is to understand conditions that lead to the existence of the
unstable 1D spiral and impact its stability in the Morris-Lecar model.
Previous studies relied on numerical simulations and shooting methods that are
impractical and inefficient at finding unstable solutions. Our computational
methods expand on techniques used to compute traveling waves and patterns as
equilibrium solutions. Therefore, after a brief description of the Morris-
Lecar system and parameters, we provide a short description of traveling waves
in reaction-diffusion equations and the computational methods for these
systems. We then end the section with a preliminary description of the
appearance of the 1D spiral.

###  2.1 Reaction-diffusion model with Morris-Lecar dynamics

Reaction-diffusion equations in one spatial dimension take the form

$\displaystyle U_{t}=DU_{xx}+F(U;\mu),\ \ U\in\mathbb{R}^{n},\
D\in\mathbb{R}^{n\times n},\ \mu\in\mathbb{R}^{p},\ x\in\mathbb{R},\
t\in\mathbb{R}$ (1)

where $U=(u_{1},u_{2},\cdots,u_{n})^{T}$ is a vector of species that diffuse
at non-negative rates given by elements of the diagonal matrix $D$. The
diffusion term couples nearby areas in space and the nonlinear reaction term
$F(U;\mu)$ defines local excitable kinematic reactions, where $\mu$ represents
the $p$ system parameters.

The Morris-Lecar model [37] is a two-variable canonical model of excitable
media, with variables $U=(V,n)^{T}$ representing the membrane potential $V$
(in mV) and the potassium channel gating variable $n$. Diffusive coupling is
only throught the membrane potential. While idealized, the model is still
framed in terms of biophysical variables and ionic currents, and parameters
can be tuned so that the it exhibits a wide range of phenomena observed in
cardiac and neural tissue. Thus, it is flexible and relatively easy to
analyze, and results can be tied to biophysical mechanisms.

Because we will be investigating the parameter impacts on the likelihood of
reflection, we give a brief overview of the physiological processes included
in the model. The voltage equation includes an inward depolarizing calcium
(Ca) current, outward rectifying potassium (K) current, and a leak ($\ell$)
current. The local reaction dynamics take the form

$\displaystyle F(V,n;\mu)=$ $\displaystyle\begin{pmatrix}f_{1}(V,n;\mu)\\\
f_{2}(V,n;\mu)\end{pmatrix}=\begin{pmatrix}-G_{Ca}m_{\infty}(V)\left(V-E_{Ca}\right)-G_{K}n\left(V-E_{K}\right)-G_{\ell}\left(V-E_{\ell}\right)\\\
\epsilon\left[\alpha(V)(1-n)-\beta(V)n\right]\end{pmatrix}$ (2)

where
$m_{\infty}(V)=0.5\left(1+\tanh\left(\frac{V-u_{1}}{u_{2}}\right)\right)$.
Values of $G_{Ca}$ and $G_{K}$ give the maximum conductances for the two
currents.

We write the reaction equations for $n$ in terms of the opening and closing
rates $\alpha(V)$ and $\beta(V)$ instead of the standard asymptotic gating
variable $n_{\infty}(V)$ and time constant $\tau_{n}(V)$ form in order to
separate the processes impacting the activation and deactivation of the
potassium channels

$\displaystyle\alpha(V)=\frac{1}{2}\left[\cosh\left(\frac{V-u_{3a}}{2u_{4a}}\right)+\tanh\left(\frac{V-u_{3a}}{2u_{4a}}\right)\cosh\left(\frac{V-u_{3a}}{2u_{4a}}\right)\right]$
$\displaystyle\beta(V)=\frac{1}{2}\left[1-\tanh\left(\frac{V-u_{3b}}{u_{4b}}\right)\right]\cosh\left(\frac{V-u_{3b}}{2u_{4b}}\right).$

To be consistent with previous studies, the excitability parameter $\epsilon$
is factored out of $\alpha$ and $\beta$. Parameters $u_{3a}$ and $u_{4b}$
shift and scale the activation rate $\alpha$, and $u_{3b}$ and $u_{4b}$ tune
the deactivation rate $\beta$. The influence of these variables on $\alpha$,
$\beta$, and the asymptotic gating function $n_{\infty}$ are show in Figure 2.
Increasing $u_{4a}$ or $u_{4b}$ corresponds to a reduction in the activation
and deactivation rates at high and low voltages, respectively. Decreasing
$u_{3a}$ shifts the activation $\alpha$ curve to lower voltages, coinciding
with the opening of the potassium gates and the activation phase of the
asymptotic variable $n_{\infty}$ at lower voltages. Variations to $u_{3b}$
have the same effect on the deactivation curve $\beta$, with increasing
$u_{3b}$ resulting in the deactivation phase of $n_{\infty}$ shifting to
higher voltages.

Figure 2: Impacts of $u_{3a}$, $u_{3b}$, $u_{4a}$, and $u_{4b}$ on (a,c)
potassium asymptotic gating variable $n_{\infty}$, and on (b,d) activation
$\alpha$ and deactivation $\beta$ rates.

Throughout this study, $\epsilon$ will be the primary bifurcation parameter,
and we also investigate the impacts of the parameters
$\mu=\\{G_{Ca},G_{K},u_{3a},u_{3b},u_{4a},u_{4b}\\}$. These parameters are
chosen as they represent changes to the key physiological processes captured
in the model, that is the depolarization and repolarization currents and the
timing of the repolarization process. Parameter values are given in the
Appendix, as well as the forms of $n_{\infty}(V)$ and $\tau_{n}(V)$.

###  2.2 Traveling Pulses and Wave Trains in Reaction-Diffusion Systems

Models of excitable media often have two traveling pulse solutions, reffered
to as ‘fast’ and ‘slow’ pulses due to their speed of propagation [25, 12].
Pulses are solutions of the form $U(x,t)=U_{f,s}(x-c_{f,s}t)$, which have a
constant profile $U_{f,s}$ that is translated at speeds $c_{f,s}$ [45, 30];
generally fast pulses are stable and slow are unstable. For rightward moving
pulses, the speeds are such that $0<c_{s}<c_{f}$ with subscripts denoting
‘fast’ or ‘slow.’ In the co-moving coordinate $\xi=x-c_{f,s}t$, the pulses are
stationary solutions of

$\displaystyle U_{t}=DU_{\xi\xi}+c_{f,s}U_{\xi}+F(U;\mu),\ \ \
\xi\in\mathbb{R}.$ (3)

In the Morris-Lecar system, pairs of pulses are created as the parameter
$\epsilon$ moves through a saddle-node bifurcation $\epsilon_{SN}$ (Figure 3).
Throughout, we will refer to $\epsilon_{SN}$ as the saddle-node bifurcation
point, and we will assume the pulses exist for $0<\epsilon<\epsilon_{SN}$.
Three basic solutions exist for values of $\epsilon$ near $\epsilon_{SN}$: the
fast and slow pulses, and the spatially homogeneous rest state. The stable
fast pulse corresponds to action potentials that propagate in the excitable
tissue and, under these conditions, the unstable slow pulse organizes the
possible outcomes and acts as a threshold solution: perturbations that bring
the system above the stable manifold of the slow pulse grow into a fast pulse,
whereas conditions that push the system below the stable manifold will decay
to the rest state.

Figure 3: (a) Bifurcation diagram in $\epsilon$ for fast and slow pulses in
Morris-Lecar. Default parameter set used. (b) Examples of fast and slow pulses
for $\epsilon=0.35$. Arrow indicates direction of motion.

Reaction-diffusion systems also support periodic traveling waves, i.e. wave
trains, which are solutions to (1) of the form $U(x,t)=U_{\infty}(\kappa
x-\omega t)$ where $U_{\infty}$ is $2\pi$-periodic in its argument and has
wave number $\kappa$ and frequency $\omega$. Wave trains arise in one-
parameter families with $\omega$ and $\kappa$ connected by the nonlinear
dispersion relation $\omega=\omega(\kappa)$ [45, 30]. In this application,
periodic wave trains make up the far-field structure of the 1D spiral wave
(see Figure 1b) and will be an important piece in the computation of the 1D
spiral wave.

The pulse and wave train solutions are accurately computed by formulating them
as equilibrium solutions in a co-moving frame [6, 20]. For example, the fast
pulse is numerically found by solving for the vector $U_{f}(\xi)$ that is a
solution to

$\displaystyle 0=DU_{\xi\xi}+c_{f}U_{\xi}+F(U;\mu),\ \ \xi\in[0,L)$

on a periodic discretized domain with analytical derivatives replaced by
fourth-order centered finite difference differentiation matrices. The speed
$c_{f}$ is a free variable and is solved for with the addition of a phase
condition that fixes the translational symmetry of the system and selects a
unique solution [10]. For pulse solutions, the domain length $L$ is taken
large enough so that the pulse does not interact with itself, and the domain
is scaled by the wave number to be $2\pi$-periodic for wave train solutions.

Stationary solutions can be numerically continued in any parameter to provide
explicit information about how parameter changes result in modifications to
the speed and pulse profile [1, 32]. Unlike direct numerical simulations,
boundary value and continuation methods can solve for both stable and unstable
patterns and saddle-node bifurcation curves. This is a crucial necessity for
our study, as the slow pulses and 1D spiral waves are unstable. Throughout, we
use secant methods to accurately trace out the nonlinear continuation curves
[10].

###  2.3 Appearance of the one-dimensional spiral wave and reflections

As $\epsilon$ is decreased below a critical value $\epsilon_{*}$, it was found
[16] that the set of solutions can become more complex as a fourth basic
solution emerges: the 1D spiral wave. The 1D spiral wave is time-periodic and
thus corresponds to an unstable periodic orbit [21, 31]; throughout we will
interchangeably refer to the pattern as an unstable periodic orbit or a 1D
spiral wave. Even though the structure is unstable, the existence of this
periodic orbit allows for new spatiotemporal dynamics, and solutions that
start near or approach the orbit may lead to a transient series of reflections
before diverging to a stable pattern. Figure 4 depicts the spatiotemporal
patterns that can arise based on the amplitude $A$ of the initial condition
and value of $\epsilon$. We emphasize that all patterns in Figure 4 are formed
on a homogenous spatial domain - no heterogeneity is inducing the reflection-
like behavior.

Figure 4: Schematic of bifurcation diagram based on amplitude $A$ of initial
condition and parameter $\epsilon$ in the Morris-Lecar model. Locations of
saddle node and heteroclinic bifurcation points indicated by $\epsilon_{SN}$
and $\epsilon_{*}$, respectively. Panels on left and right show resulting
spatiotemporal pattern for indicated position along left and right vertical
lines. Simulations run on a homogeneous domain.

The branch of 1D spiral solutions appears to bifurcate from the branch of slow
pulse solutions as illustrated in Figure 4, and it is has been found that the
period of the 1D spiral increased as $\epsilon$ approached $\epsilon_{*}$.
This led Cytrynbaum & Lewis to postulate that the 1D spiral arises at
$\epsilon_{*}$ though a rearrangement of heteroclinic connections of the
unstable slow pulse in a global bifurcation [16].

Our goals are to determine what conditions control the existence of the
unstable periodic orbit underlying the 1D spiral pattern and factors that
encourage and depress reflections when the 1D spiral does exist. The existence
of 1D spiral, including the bifurcation structure giving rise it, is
investigated in Section 4, but we first address the likelihood of reflections
in Section 3.

##  3 Parameter-dependent propensity for reflections

First, we consider scenarios in which 1D spirals exist
($\epsilon<\epsilon_{*}$) and measure the parameter-dependent propensity for
reflections. Specifically, we compute the 1D spiral as a spatiotemporal
equilibrium solution and consider its spectral stability, as the magnitude of
the unstable eigenvalue informs how repelling the periodic orbit is. The more
unstable the periodic orbit, the closer a system must approach the 1D spiral’s
stable manifold to produce a reflection. In this section, we use direct
numerical simulations to link conditions that produce reflections to the
magnitude of the unstable eigenvalue. Then, we investigate parameter-
dependence of the eigenvalue and determine conditions that promote or deter
reflection behavior.

###  3.1 Reflections induced by elevated refractory region

Reflections can be generated by a stable fast pulse interacting with a
functional heterogeneity, such as a region with a transiently elevated
refractory level. To demonstrate that altering system conditions influences
the likelihood of reflections and to motivate the need for more precise
computational methods, we test the system’s propensity for reflections using
direct numerical simulation.

The numerical experiments are performed as follows. A fast pulse is elicited
so that it travels toward a small region with an elevated refractory level.
First, the pulse is evolved for $t=5$ to ensure that a stable traveling fast
pulse solution is reached. A refractory bump of the form

$\displaystyle B\exp\left(\frac{-(x-x_{0})^{2}}{\sigma^{2}}\right)$

is added to the recovery variable $n$, where $x_{0}$ is set to 0.15 in front
of the recovery peak and $\sigma=0.05$. The solution is then evolved for
$t=10$, to ensure that the pulse interacts with the recovery bump. For
sufficiently small amplitude $B$, refractory bumps cause the fast pulse to
slow slightly, but the pulse ultimately propagates through the region and
returns to the original stable form of the fast pulse. For large amplitude
$B$, the pulse propagation is blocked and annihilated. Amplitudes between
these two cases can generate one or more reflections (Figure 5). We define
$[B_{\text{min}},B_{\text{max}}]$ as the interval of bump amplitudes that led
to at least one reflection. The interval endpoints are computed via a
bisection method outlined in the Appendix. In cases where reflections did
occur, the length of the bump amplitude interval
$B_{\text{max}}-B_{\text{min}}$ was taken to be one measure of the propensity
for reflections for the given parameter set.

Figure 5: (a) Fast pulse interacting with a refractory bump, indicated with
green arrow. (b) Series of space-time plots showing the voltage variable of a
stable false pulse interacting with a refractory bump of increasing amplitude.
From left to right amplitudes $B$ are 0.349, 0.34957, and 0.35.

Using the experimental design described above, we test how varying the
potassium conductance $G_{K}$ and shifting the potassium deactivation process
$u_{3b}$ alters the likelihood of reflections. We find that increasing
$G_{K}$, which corresponds to strengthening the repolarization current,
results in a decreased tendency for reflections (Figure 6). Likewise, the bump
amplitude interval and therefore propensity for reflections decreases with
$u_{3b}$.

Figure 6:  Interval of refractory bump amplitudes that lead to reflection as a
function of parameter (a) $G_{K}$ and (b) $u_{3b}$.

###  3.2 Computation of 1D spiral waves

The results of the refractory bump experiment yields information for how one
form of heterogeneity produces reflections by perturbing the stable pulse so
that the system is near the unstable 1D spiral. To efficiently determine the
propensity of reflections across all types of heterogeneities and various
parameters, we directly compute the 1D spiral and investigate its spectra.
Specifically, we focus on the magnitude of the unstable eigenvalue because it
dictates how repelling the unstable periodic orbit is; systems with a larger
positive eigenvalue will diverge from the periodic orbit quicker and only
initial conditions very close to the orbit will stay near it long enough to
generate a transient series of reflections. Therefore, it is easier to
generate reflections when the 1D spiral is less unstable (lower value of the
unstable eigenvalue). To compute the 1D spiral and its stability, we consider
it as a spatiotemporal pattern and explicitly take advantage of the pattern’s
form.

The 1D spiral has the structure of an anti-symmetric spatiotemporal source
defect [46]: the core is a source that emits a series of counter-propagating
fast pulses that assemble in the far-field into a regular structure of
periodic wave trains with wave number $\kappa$ and frequency $\omega$ uniquely
selected by the core. The 1D spiral waves are time $T$-periodic solutions
$U_{*}(x,t)$ of equation (1) with $U_{*}(x,t+T)=U_{*}(x,t)$.

On the two-dimensional spatiotemporal domain $\Omega=[-L,L]\times S^{1}$, the
1D spiral wave can be considered as an equilibrium solution of

$\displaystyle 0=DU_{xx}-\omega U_{\tau}+F(U;\mu),\ \ (x,\tau)\in[-L,L]\times
S^{1}$ (4)

where $\tau=\omega t$, $\omega=2\pi/T$, and boundary conditions at $x=\pm L$
are chosen to mimic an infinite spatial domain. By utilizing the core and far-
field structure of the pattern, the 1D spiral wave can be written as

$\displaystyle U_{*}(x,\tau)=W(x,\tau)+\chi(x)U_{\infty}(\kappa x-\tau)$

where the smooth cut-off function
$\chi(x)=1-0.5\left(\tanh[m(x+d)]-\tanh[m(x-d)]\right)$ connects the far-field
wave train $U_{\infty}(\kappa x-\tau)$ with the core solution $W(x,\tau)$ [34,
22, 19]. The far-field solution is pre-computed by solving for the wave train
on a 1D periodic domain and interpolating onto the 2D spatiotemporal domain
$\Omega$. Thus, when solving the system, the root-finding framework exploits
the core and far-field structure of the pattern for added efficiency and only
solves for the localized (core) function $W(x,\tau)$.

Because the 1D spiral is formulated as an equilibrium solution, its stability
can be determined from the spectra of the operator $\mathcal{L}_{*}$ formed by
linearizing (4) about the 1D spiral $U_{*}(x,\tau)$

$\displaystyle\mathcal{L}_{*}V=DV_{xx}-\omega
V_{\tau}+F_{U}(U_{*}(x,\tau);\mu)$

and considering the resulting eigenvalue problem $\mathcal{L}_{*}V=\lambda V$
for $\lambda\in\mathbb{C}$.

Formulated in the above way, computing the 1D spiral as a spatiotemporal
pattern corresponds to directly finding the unstable periodic orbit.
Furthermore, we locate a single, real unstable eigenvalue $\lambda_{u}$ with
corresponding eigenfunction localized near the core (shown in Appendix); the
value of this unstable eigenvalue informs how repelling the orbit is. We claim
reflections will be harder to generate when the eigenvalues are more positive
(i.e., the system is more unstable). To support this claim, we repeat the
refractory experiment for all parameters $\mu$ of interest and compare the
length of the amplitude intervals that generate reflections
($B_{max}-B_{min}$) to the unstable eigenvalues of the 1D spiral. Figure 7
shows that the two quantities are correlated, corroborating the hypothesis
that the 1D spiral is the fundamental dynamical mechanism underlying the
ability for pulses to reflect. Thus, the unstable eigenvalue $\lambda_{u}$ of
the 1D spiral can be used as a metric for characterizing the relative
likelihood of reflections created by any heterogeneity, with more positive
eigenvalues corresponding to a lower chance of reflection.

Figure 7:  Interval of refractory bump amplitudes $B$ that lead to reflections
versus the unstable eigenvalue $\lambda_{u}$ of the 1D spiral wave. Note the
log-scale of the vertical axis.

###  3.3 Parameter dependent stability of 1D spiral and likelihood of
reflections

Under the assumption that the periodic orbit exists, we seek to understand
what system conditions contribute to a higher occurrence of reflection
behavior. To do so, the unstable eigenvalue $\lambda_{u}$ is continued in
parameters $\mu$ for fixed $\epsilon=0.2$. Continuation results are reported
in Figure 8. To aid in cross-parameter comparison, the eigenvalue continuation
curves are shown as a function of percent change in each parameter, with the
axis capped at 100%. The unstable eigenvalue is monotonic in parameters
$G_{K}$, $u_{3a}$, $u_{3b}$, and $u_{4b}$, and the periodic orbit becomes more
unstable with decreasing $u_{3a}$ and $u_{3b}$ and increasing $G_{K}$ and
$u_{4b}$. The eigenvalue behavior with respect to $G_{Ca}$ and $u_{4a}$ is
non-monotonic, with both increases and decreases in these parameters further
destabilizing the 1D spiral.

Figure 8: Results from continuation of the unstable eigenvalue $\lambda_{u}$
in each parameter. Displayed as a function of percent change in each
parameter. Fixed $\epsilon=0.2$.

We also measure the sensitivity of the unstable eigenvalues to each parameter
by (1) determining the percent change in $\lambda_{u}$ when $\mu$ is increased
by 10%, and (2) computing the maximum partial derivatives of the eigenvalue
continuation curves with respect to the corresponding parameter over each
parameter’s tested range. The percent change indicates the eigenvalue
sensitivity to small changes near the default values and the maximal partial
derivative provides a sensitivity index across the full tested parameter
range. Both values are reported in Table 1; negative values in the partial
derivative indicate the maximum value occurs for a decrease in the parameter
value. The overall highest reduction of reflection likelihood comes from
strengthening the potassium current by increasing $G_{K}$
($\max\frac{\partial\lambda_{u}}{\partial G_{K}}=8.06$). Combined, the two
sensitivity indices indicate that the stability of the periodic orbit is most
sensitive to variations in the conductances $G_{Ca}$ and $G_{K}$, moderately
sensitive to $u_{3a}$ and $u_{4b}$, and relatively insensitive to $u_{3b}$ and
$u_{4a}$. Consequences of these findings for reducing reflections and links to
the physiological processes will be further discussed in Section 6.

Table 1: Parameter dependent sensitivity of the unstable eigenvalue $\lambda_{u}$. Column 2 gives percent change in $\lambda_{u}$ corresponding to a 10% change in $\mu$. Column 3 reports the maximum value of the partial derivative of $\lambda_{u}$ continuation curve over each parameter’s tested range.  Parameter $\mu$ | % Change $\lambda_{u}$ | $\max\frac{\partial\lambda_{u}}{\partial\mu}$
---|---|---
$G_{Ca}$ | 5.69 | -2.04
$G_{K}$ | 5.46 | 8.06
$u_{3a}$ | -0.71 | -0.29
$u_{3b}$ | -0.13 | -0.07
$u_{4a}$ | -0.28 | -0.03
$u_{4b}$ | 1.42 | 0.1

##  4 Parameter-dependence for existence of the 1D spiral

Here, we investigate the bifurcation that leads to the formation of the 1D
spiral wave and how the 1D spiral’s existence depends on system parameters.
Our focus is on uncovering mechanisms that limit the 1D spiral’s existence in
order to better understand conditions that hinder reflections. We start by
confirming that the unstable periodic orbit does not originate via a local
bifurcation of the slow pulse and follow this with a description of the
proposed global bifurcation through which it arises. We then exploit the
heteroclinic structure and develop a computation method to directly compute
the global bifurcation point. Finally, we explore the parameter-dependence for
the existence of the 1D spiral and conditions that lead to the formation of
the 1D spiral are identified.

###  4.1 Evidence of a global bifurcation

The 1D spiral waves are hypothesized to arise from a global rearrangement of
heteroclinic connections [16]. To verify that the unstable periodic orbit
underlying the 1D spiral does not originate through a local bifurcation from
the slow pulse, we check spectrum of the slow pulse. The spectrum of a pulse
consists of discrete point eigenvalues and curves of continuous spectrum;
instabilities in either of these sets results in a local bifurcation. Near the
saddle-node bifurcation at $\epsilon_{SN}$, the slow pulse has one real
unstable point eigenvalue and a one-dimensional unstable manifold. We compute
the discrete spectrum of the slow pulse to determine if any additional
unstable eigenvalues appear. For $\epsilon<\epsilon_{SN}$, fast and slow
pulses were computed as solutions of (3) and numerical continuation was used
to trace out the saddle node bifurcation in $\epsilon$ (Figure 3). Eigenvalues
in the discrete spectrum of the slow pulse were computed by linearizing the
reaction-diffusion operator about the slow pulse [6]. We find that spectrum of
the slow pulse only contains one unstable point eigenvalue for
$0<\epsilon<\epsilon_{SN}$; the eigenvalue corresponding to the saddle-node
bifurcation.

Instabilities of the continuous spectrum are associated with instabilities of
the asymptotic rest state, in this case corresponding to the resting
potential. Source defects have been shown in some instances to be the result
of a continuous instability [46], however, this is not the case in the
situation presented here. The essential spectrum was computed using
continuation methods described in [42] and we find that the essential spectrum
of the slow pulse remains stable for all $0<\epsilon<\epsilon_{SN}$
(eigenvalue computations in the Appendix). Because the slow pulse does not
have any spectral changes, the 1D spiral cannot form as a spectral bifurcation
of the slow pulse.

###  4.2 Proposed global bifurcation structure

To introduce the mathematical setting and proposed heteroclinic bifurcation,
we first summarize the behavior of a set of counter-propagating of pulses on
the real line for values of $\epsilon$ near the bifurcation. Following the
notation in [16], we denote a series of traveling pulses using $f,s$ to
signify leftward fast and slow pulses, and $F,S$ to signify rightward
propagating fast and slow pulses. The schematic in Figure 9 summarizes the
proposed heteroclinic bifurcation. For $\epsilon>\epsilon_{*}$, when
perturbed, slow pulses ($S$) either decay to the rest state ($R$) or grow into
a fast pulse ($F$); that is, the unstable manifold of the slow pulse
$W^{u}(S)$ connects to $R$ and $F$. Similar behavior is seen for a slow pulse
that lies in a series of pulses. For example, an initial condition near $fSF$
(left fast, right slow, right fast) will evolve into either $fF$ or $fFF$;
that is, the unstable manifold of $fSF$ connects to $fF$ and $fFF$. As
$\epsilon$ passes through the bifurcation point $\epsilon_{*}$, it is
hypothesized that the unstable manifold of the slow pulse switches from being
connected to the fast pulse $F$ to a pair of counter-propagating fast pulses
($fF$). Now, for initial conditions near $fSF$, the system will evolve into
$fF$ if the slow pulse decays to rest, or evolve into $ffFF$ if the slow pulse
splits into counter-propagating fast pulses. It is this rearrangement of
unstable manifolds that gives rise to the unstable periodic orbit (thick black
line in Figure 9). That is, the 1D spiral’s formation relies on the
simultaneous rearrangement of all heteroclinic connections that are linked to
$W^{u}(S)$.

Figure 9: Schematic of proposed heteroclinic bifurcation structure as
$\epsilon$ passes through $\epsilon_{*}$. Throughout, empty triangles
represent unstable solutions; solid circles and squares are stable. For
$\epsilon\gg\epsilon_{*}$, unstable manifolds of the slow pulse (red arrows)
connect to the fast pulse $F$ and rest state $R$. At $\epsilon=\epsilon_{*}$,
unstable manifold of the slow pulse connects to a manifold of an oppositely
propagating slow-fast pulse pair (sF). Unstable periodic orbit (1D spiral)
appears for $\epsilon<\epsilon_{*}$ (thick black curve) as unstable manifold
of slow pulse connects to an oppositely propagating pair of fast pulses.
Dashed lines show one-dimensional stable manifold of the slow pulse. Figure
reproduced with minor revisions and permission from [16].

During this rearrangement of unstable manifolds, the temporal period of the 1D
spiral will limit to infinity as $\epsilon\rightarrow\epsilon_{*}$, and at the
critical point $\epsilon_{*}$, the limiting behavior along $W^{u}(S)$ is a
rightward propagating slow pulse ($S$) that splits into a rightward fast pulse
and leftward slow pulse ($sF$). This behavior corresponds to a heteroclinic
connection from $S$ to $sF$, or equivalently the unstable manifold of the slow
pulse $W^{u}(S)$ to a manifold of an counter-propagating slow-fast pulse pair
$M_{sF}(d)$, with solutions on the manifold parameterized by the distance $d$
between the fast-slow pulse peaks. (We will exploit this structure to compute
the global bifurcation point $\epsilon_{*}$ in sections 4.4 and 4.5.)

###  4.3 Further evidence for a heteroclinic bifurcation

To support that the global bifurcation is of the heteroclinic type, we
consider the behavior of the temporal period of the 1D spiral as $\epsilon$
increases to the bifurcation point $\epsilon_{*}$. Using the methods described
in Section 3.2, 1D spiral waves are formulated as a spatiotemporal equilibrium
pattern on the domain $\Omega=[-1,1]\times S^{1}$ and numerically computed.
Continuation of the pattern in parameter $\epsilon$ provides information on
how the temporal period $T$ and far-field spatial wave number $\kappa$ change
with parameters (Figure 10). The computations illustrate that the temporal
period $T$ scales like $\log\left(\epsilon_{*}-\epsilon\right)$, which is the
expected scaling of a heteroclinic bifurcation [47]. Consequently, as
$\epsilon\rightarrow\epsilon_{*}$ we find that the perturbed slow pulse at the
core of a 1D spiral persists for a longer time before destabilizing into a
fast pulse and retrograde slow pulse, which is consistent with the proposed
bifurcation structure described above.

Other system parameters $\mu$ can be used in place of $\epsilon$ when
considering the 1D spiral existence interval and bifurcation point.
Continuations of the 1D spiral in the parameters
$\mu=\\{G_{Ca},G_{K},u_{3a},u_{3b},u_{4a},u_{4b}\\}$ show that the periodic
orbit terminates at a point $\mu_{*}$, with the period $T$ scaling like
$\log(\mu_{*}-\mu)$ for parameter values $\mu$ near $\mu_{*}$ (Figure 10).
These results suggest the heteroclinic bifurcation that results in the 1D
spiral is a robust and more general phenomenon than just observed here.

Figure 10:  Morris-Lecar Model. (a) Scaling of temporal period $T$ as solution
is numerically continued in $\epsilon$ to $\epsilon_{*}$. (b) Spatiotemporal
patterns corresponding to indicated positions on the continuation. Time $t$
rescaled to laboratory frame. (c) Temporal period $T$ versus
$\log\left(\mu-\mu_{*}\right)$ for parameters $\mu$ near $\mu_{*}$.

###  4.4 Computation of the heteroclinic bifurcation point

We develop a numerical method to directly compute the heteroclinic connection
between the slow pulse ($S$) and the fast-slow counter-propagating solution
($sF$), in order to compute the global bifurcation point $\epsilon_{*}$ that
gives rise to the 1D spiral. We provide a summary and justification for the
method here, and describe implementation details in the appendix. Heteroclinic
connections can be numerically computed as a truncated boundary value problem
on a finite time domain by augmenting the system with projection boundary
conditions that ensure the orbit approaches the equilibria in finite time
[10]. We extend these methods to directly compute the heteroclinic connection
between the unstable manifold of the slow pulse $W^{u}(U_{s})$ and manifold of
an oppositely propagating slow-fast pair $M_{sF}(d)$. Figure 11 shows a
schematic of the proposed heteroclinic orbit at $\epsilon_{*}$ and
corresponding spatiotemporal pattern that connects $S$ to $M_{sF}(d)$.
Specifically, we formulate the heteroclinic connection as an equilibrium
solution on a spatiotemporal domain in the co-moving frame of the slow pulse

$\displaystyle 0=DU_{\xi\xi}+c_{s}U_{\xi}+F(U;\mu)-U_{t},\ \
(\xi,t)\in[-L,L]\times[-T,T]$

equipped with the projected boundary conditions

$\displaystyle U$ $\displaystyle(\xi,t=-T)=U_{s}(\xi)+\gamma V_{s}^{u}(\xi)$
(5) $\displaystyle 0$ $\displaystyle=\langle
M(\xi)-U(\xi,t=T),\psi^{*}(\xi)\rangle.$ (6)

Condition (5) implies that the solution at $t=-T$ is the slow pulse
$U_{s}(\xi)$ plus a small perturbation $\gamma$ in the direction of the
unstable eigenfunction $V_{s}^{u}(\xi)$. Condition (6) enforces that the
difference between the solution at $t=T$ and a solution $M(\xi)$ on the
manifold $M_{sF}(d)$ is perpendicular to the unstable adjoint eigenfunction
$\psi^{*}(\xi)$ of $M(\xi)$. The necessary pulses, eigenfunctions, and adjoint
eigenfunctions are additionally solved as equilibrium solutions on 1D domains,
and secant continuation allows the bifurcation point $\epsilon_{*}$ to be
continued in system parameters.

In defining this method, we assume that the unstable manifold of $M_{sF}$ is
one-dimensional, and more specifically that $W^{U}(M_{sF})$ corresponds to
$W^{u}(U_{s})$. This assumption is supported by results in [51] that prove
systems with traveling pulse solutions also support solutions that are nearly
the linear supposition of oppositely propagating pulses. Under this
formulation, the eigenspaces of the counterpropagating pulses are found to be
a cross-product of the eigenspaces of each individual pulse. Because
$W^{U}(M_{sF})$ and $W^{u}(U_{s})$ are both one-dimensional, heteroclinic
connections from $U_{S}$ to $M_{sF}$ only occur at isolated parameter values
[26]; for each set of parameters $\mu$, this system has one unique solution
that corresponds to the global bifurcation point $\epsilon_{*}$. We include
$\epsilon_{*}$ as an unknown parameter in the computations and directly solved
for its value.

Figure 11: (a) Schematic of the hypothesized heteroclinic connection between
$U_{s}$ and $M_{sF}$ at $\epsilon=\epsilon_{*}$. Truncated boundary value
problem connects a solution that starts near $U_{s}$ to a solution that ends
near $M_{sF}$. (b) Form of spatiotemporal pattern associated with the
heteroclinic connection, in co-moving frame of initial slow pulse.

###  4.5 Parameter-dependent sensitivity of the Bifurcation Points
$\epsilon_{*}$ and $\epsilon_{SN}$

A major goal in this study is to assess the parameter dependent sensitivities
of the heteroclinic bifurcation point $\epsilon_{*}$, since shifting
$\epsilon_{*}$ will alter the existence window for reflections. However,
changes to system parameters will also impact $\epsilon_{SN}$ and the
existence of traveling pulses, thus we monitor the values of both
$\epsilon_{*}$ and $\epsilon_{SN}$. Continuations of the points $\epsilon_{*}$
and $\epsilon_{SN}$ in $G_{Ca}$, $G_{K}$, $u_{3a}$, $u_{3b}$, $u_{4a}$, and
$u_{4b}$ establish the independent influences of each parameter on the
bifurcation points. Values of $\epsilon_{*}$ and $\epsilon_{SN}$ were
continued in each parameter until, at minimum, either $\epsilon_{*}$ fell
outside the range $[0.2,0.4]$ or the continuation parameter changed by at
least $\pm$100%.

We find that $\epsilon_{*}$ and $\epsilon_{SN}$ change monotonically with all
parameters except $u_{4a}$. Specifically, increasing $G_{Ca}$, $u_{3a}$, and
$u_{3b}$ increase both bifurcation points, and increasing $G_{K}$ and $u_{4b}$
decreases $\epsilon_{*}$ and $\epsilon_{SN}$. Behaviors observed with changes
in $u_{4a}$ are more complex. Small increases in $u_{4a}$ from the default
value lead to decreases in $\epsilon_{*}$ and $\epsilon_{SN}$, but larger
increases eventually decrease both $\epsilon_{*}$ and $\epsilon_{SN}$; these
behaviors will be further discussed below.

The sensitivity of the bifurcation points is quantified similarly to that of
$\lambda_{u}$ in Section 3 by computing: (1) the percent change in
$\epsilon_{*}$ and $\epsilon_{SN}$ corresponding to 10% increases in parameter
values, and (2) the maximum partial derivative of $\epsilon_{*}$ and
$\epsilon_{SN}$ continuations with respect to each parameter. These values
present two sensitivity measures of the bifurcation points to local changes in
the parameter values and are listed in Table 2. From these two values, we find
that $\epsilon_{*}$ is most sensitive to changes in the calcium conductance
$G_{Ca}$; this sensitivity is an order of magnitude higher than the other
parameters in terms of the partial derivative. The potassium conductance
$G_{K}$ and $u_{3a}$ have moderate influence on $\epsilon_{*}$, followed by
$u_{4b}$, $u_{4a}$, and $u_{3b}$.

Table 2: Parameter-dependent changes in $\epsilon_{*}$ and $\epsilon_{SN}$. Columns 2 and 3 are the percent change after a 10% increase in parameter. Columns 4 and 5 are the maximum value of the partial derivative of each bifurcation point continuation over each parameter’s tested range. Negative values indicate increasing the parameter decreases the bifurcation location.  Parameter $\mu$ |  $\%$ Change $\epsilon_{*}$ |  $\%$ Change $\epsilon_{SN}$ | $\max\frac{\partial\epsilon_{*}}{\partial_{\mu}}$ | $\max\frac{\partial\epsilon_{SN}}{\partial_{\mu}}$
---|---|---|---|---
$G_{Ca}$ | 94.36 | 63.18 | 0.827 | 0.722
$G_{K}$ | -13.9 | -10.82 | -0.089 | -0.097
$u_{3a}$ | 3.04 | 2.26 | 0.089 | 0.069
$u_{3b}$ | 0.50 | 0.64 | 0.018 | 0.011
$u_{4a}$ | -4.32 | 1.15 | -0.032 | 0.018
$u_{4b}$ | -3.35 | -0.18 | -0.058 | -0.028

The ratio $\epsilon_{*}/\epsilon_{SN}$ provides a measure of the existence
interval of the 1D spiral relative to that of the traveling pulse. Ratios near
one represent conditions where reflections are possible for a large portion of
the pulse existence window. The relative bifurcation values were tracked for
the same changes in parameters (Figure 12) and the ratio
$\epsilon_{*}/\epsilon_{SN}$ is monotonic in all parameters. Decreases in
$G_{Ca}$, $u_{3a}$, and $u_{3b}$ and increases $G_{k}$, $u_{4a}$, and $u_{4b}$
lead to reductions in $\epsilon_{*}/\epsilon_{SN}$. The most significant
decreases correspond to increasing $u_{4a}$ and $u_{4b}$.

The ratio $\epsilon_{*}/\epsilon_{SN}$ does not fully characterize the
relative length of the 1D spiral existence interval as both bifurcation points
shift. It is also important that $\epsilon_{SN}$ remains high so that
traveling pulses exist over a large range of maximal excitabilities.
Corresponding values of $\epsilon_{*}$ and $\epsilon_{SN}$ are summarized in
Figure 12b. For example, while variations in $G_{Ca}$, $G_{K}$, and $u_{3a}$
reduce $\epsilon_{*}/\epsilon_{SN}$, these reductions in the reflection
existence interval coincide with a similar reduction in the traveling pulse
existence interval. These findings are consistent when comparing the
sensitivities of $\epsilon_{*}$ and $\epsilon_{SN}$ in Table 2. The saddle
node bifurcation point $\epsilon_{SN}$ has a similar level of sensitivity as
$\epsilon_{*}$ to the conductances $G_{Ca}$ and $G_{K}$ and also $u_{3a}$.

More favorable outcomes (i.e., those with a decreasing ratio and relatively
high $\epsilon_{SN}$) are accomplished by increases to $u_{3b}$, $u_{4a}$, and
$u_{4b}$. Likewise, the sensitivity index of $\epsilon_{*}$ to the scaling
parameters $u_{4a}$, and $u_{4b}$ is twice as great as $\epsilon_{SN}$.
Moderate increases to $u_{4a}$ precipitate the most advantageous scenario: a
decreasing 1D spiral existence interval coupled with a growing
$\epsilon_{SN}$. Connections to the physiological processes will be revisited
and discussed in Section 6.

Figure 12: (a) Ratio $\epsilon_{*}/\epsilon_{SN}$ as function of parameter
variations. (b) Relative locations of $\epsilon_{*}$ and $\epsilon_{SN}$.

##  5 Type I vs type II excitable systems

In 1996, Ermentrout and Rinzel [21] used one-dimensional tissue models and
systems of coupled neurons to investigate reflection behaviors. They showed
that the existence of an unstable periodic orbit in a phase model of two
coupled cells whose dynamics mimicked a system with a saddle-type threshold,
i.e., type I excitability. From their results, they postulated that, while
reflections may be possible in systems with type II excitability, they more
robustly occur over a wider range of conditions in type I systems. Later
results reported that reflections could not occur in type I systems under some
conditions, specifically those near the traveling pulse saddle-node
bifurcation [16], and 1D spirals have been detected in a system with type II
FitzHugh-Nagumo dynamics [15]. To investigate Ermentrout and Rinzel’s proposed
theory for the influence of type I and type II local dynamics on the
susceptibility to reflections, we interpret our results on the 1D spiral’s
existence and stability from Sections 3 and 4 in the context of the
excitability type. We are particularly interested in understanding if and how
type I and type II excitability influences the tendency of reflections.

An excitable system can be classified as type I or II excitable depending on
the mechanisms facilitating the transition from excitable to oscillatory
dynamics in the local reaction terms and the structure of the nullclines (see
Figure 13). Type I systems have three equilibrium points, which we denote $A$,
$B$, and $C$ in order of increasing voltage. Point $A$ is a stable node, $B$
is a saddle point, and $C$ is an unstable spiral or node. The stable manifold
of $B$ acts as a strict firing threshold. Ermentrout and Rinzel postulated
that it is this true threshold and subsequent consistent large amplitude
action potentials characteristic of type I systems that promote robust
reflection phenomena [21]. The transition from type I to type II excitability
is defined by the collision and annihilation of the unstable depolarized
equilibria $B$ and $C$ via a saddle-node bifurcation. Type II systems thus
have $A$ as the single equilibrium point and lack the saddle structure; these
excitable systems instead display a quasi-threshold behavior.

Figure 13: Nullclines in the Morris-Lecar model demonstrating (a) Type I
excitability and (b) Type II excitability. Equilibrium points labeled with
$A$, $B$, and $C$.

To differentiate amongst type I conditions, we classify systems with a larger
Euclidean distance between the depolarized equilibria $B$ and $C$ as more
“type I-like”. Likewise, type II excitable systems will be quantified by the
minimum distance between the $n$\- and $V$\- nullclines, restricted to
voltages between -30 and 0, i.e., near the $B$-$C$ saddle-node bifurcation.
Conditions that lead to a larger separation in the nullcines will be
considered as more “type-II like”. Figures 14c-d demonstrate the transition
from type I to type II excitability in a parameter that modifies the
$n$-nullcline ($u_{4a}$) and one that alters the $V$-nullcline ($G_{K}$).

Figure 14a plots the unstable eigenvalues as a function of the type I and type
II distances. Values on the horizontal-axis represent the level of type I to
type II-ness of the systems, with positive values indicating the distance
between equilibria $A$ and $B$ (type I) and negative values are the distance
between the $n$ and $V$ nullclines (type II). Note that with the default
parameter values, the system is of type I. Eigenvalue continuations were
performed with $\epsilon=0.1$ for $G_{Ca}$ and $G_{K}$ to allow the 1D spiral
to persist into the type II parameter range. All other eigenvalue
continuations kept $\epsilon=0.2$. Type II excitability was not achieved
through significant variations in parameter $u_{4b}$ since $u_{4b}$ modifies
the top bend of the $n$-nullcline with only minimal impact on the bottom bend.
Likewise, due to the $\epsilon$-existence interval of the 1D spiral,
eigenvalue continuations of $G_{Ca}$ and $u_{3a}$ only minimally crossed into
the type II region.

Overall, the trends indicate when conditions are more type I-like, the
unstable periodic orbit is less repelling, hence reflections are more likely
to occur in the system. The opposite is true for type-II like systems; the
further the separation in the $n$ and $V$ nullclines, the more repelling the
periodic orbit is, which would result in a lower probability of reflections.
Figure 8 shows that parameters $G_{Ca}$ and $u_{4a}$ exhibit non-monotonic
trends in the eigenvalues, yet we do find that type II local dynamics
generally result in a more unstable 1D spiral. Across all parameters, the
$\epsilon_{*}/\epsilon_{SN}$ ratio monotonically decreases as the system
transitions from type I to type II excitable dynamics, again indicating
reflections are less likely as they exist for a smaller portion of the
traveling wave interval.

Our findings support Ermentrout and Rinzel’s claim that local dynamics of type
I excitability are more prone to reflections, however the results are not
clear-cut. Type I systems tend to have a less repelling unstable periodic
orbit, but as the eigenvalue continuations show (Figure 8) there are cases
where unstable eigenvalue response is not monotonic in the shift from type I
to type II. Furthermore, 1D spirals tend to exist for a larger portion of the
traveling pulse existence interval in type I systems, but the decreasing
$\epsilon_{*}/\epsilon_{SN}$ ratio often coincides with a reduced traveling
wave existence interval. We furthermore find that 1D spirals exist over a
broad range of parameters in type II excitable systems, and therefore
reflections are possible.

Figure 14: (a) Unstable eigenvalue continuations as a function of the type I
and type II distances. Positive values on the horizontal axis indicate a
larger distance between the two depolarized equilibria (b)
$\epsilon_{*}/\epsilon_{SN}$ ratios as a function of type I and type II
excitability. (c) Phase planes with varying values of $u_{4a}$. (d) Phase
planes with varying values of $G_{K}$.

##  6 Discussion

Life-threatening arrhythmia can arise when action potentials generate
reflections that lead to a localized backward propagating pulse, thus it is
crucial to elucidate the mathematical structure underlying reflection
behavior. We utilized the Morris-Lecar model to demonstrate that
distinguishing conditions that increase the propensity for generating
reflections is equivalent to identifying conditions that promote the existence
and stability of the unstable periodic orbit corresponding to the 1D spiral
wave.

In the Morris-Lecar system, fast and slow traveling pulses exist for values of
$\epsilon$ below the saddle-node bifurcation at $\epsilon_{SN}$. For
$\epsilon$ near $\epsilon_{SN}$, the slow pulse acts as the threshold
separating the basins of attraction for the stable fast pulse and rest state.
The possible wave behaviors become more complex for
$\epsilon<\epsilon_{*}<\epsilon_{SN}$ as the 1D spiral wave appears. Even
though the 1D spiral is an unstable structure, its existence is significant
because it introduces the possibility of reflections.

The 1D spiral wave was proposed to emerge through a rearrangement of the
unstable manifold of the slow pulse in a heteroclinic bifurcation at
$\epsilon_{*}$ [16]. Our findings support the proposed heteroclinic
bifurcation structure. First, we confirmed the slow pulse does not undergo any
spectral changes for $0<\epsilon<\epsilon_{SN}$, and additionally found that
the 1D spiral period scales like $\log(|\epsilon-\epsilon_{*}|)$, as expected
for a heteroclinic bifurcation. The central element of the heteroclinic
bifurcation is that at $\epsilon=\epsilon_{*}$, the unstable manifold of the
slow pulse rearranges and connects to a manifold of a counter-propagating
slow-fast pulse pair. We directly computed this connection numerically as a
truncated boundary value problem to solve for the value of $\epsilon_{*}$ and
determined how $\epsilon_{*}$ changed with conditions. The computation was
crucial in determining the parameter-dependence of $\epsilon_{*}$, but it also
justifies the proposed structure.

Numerical computations additionally provided insights into the stability of
the 1D spiral. Framing the 1D spiral as a spatiotemporal equilibrium solution
allowed us to compute the spectra of the linearized operator. The single
positive eigenvalue of the 1D spiral informed how repelling the unstable
periodic orbit was, and yielded a measure of reflection propensity.

To minimize the possibility of reflections in the Morris-Lecar model, the
optimal physiological scenario is a stable fast pulse solution existing for a
large range of $\epsilon$ and a highly unstable 1D spiral existing only on a
narrow interval. This optimal scenario corresponds to conditions that have a
high value of $\epsilon_{SN}$, a low $\epsilon_{*}/\epsilon_{SN}$ ratio, and a
highly unstable eigenvalue for $\epsilon<\epsilon_{*}$. Even if optimal
conditions are not likely to be fully achieved, it may be possible to suppress
reflections by either shifting the system outside of the 1D spiral existence
interval, or further destabilizing the periodic orbit.

For conditions in which the 1D spiral exists, the parameters $G_{K}$,
$G_{Ca}$, $u_{4b}$ and $u_{3a}$ have the most influence on the unstable
eigenvalue of the 1D spiral. The magnitude of the unstable eigenvalue of the
1D spiral captures how repelling the unstable solution is, and therefore how
close the system must be to the 1D spiral in order to generate a reflection.
Increasing $G_{K}$ and $u_{4b}$ leads to the largest increases in the positive
eigenvalue, whereas decreases in $u_{3a}$ and $G_{Ca}$ lead to more moderate
eigenvalue increases. These parameter variations further destabilize the
periodic orbit and reduce the likelihood of reflections as only a very narrow
range of conditions will bring the system close enough to the periodic orbit
to produce reflections.

Changes in $u_{3a}$, $u_{4a}$, and $u_{4b}$ most significantly decrease
$\epsilon_{*}$ while maintaining a large existence window of traveling pulses.
Specifically, increasing $u_{4a}$ or decreasing $u_{3a}$ show large decreases
in $\epsilon_{*}$ relative to $\epsilon_{SN}$. Modifications to the maximum
calcium and potassium conductances ($G_{Ca}$ and $G_{K}$) greatly impact
$\epsilon_{*}$, but also lead to similar transformations in $\epsilon_{SN}$.
In fact, in these cases it is unclear if the parameter changes decrease
$\epsilon_{*}$ directly, or if it is simply a byproduct of the decrease in
$\epsilon_{SN}$. Thus, while these parameters have substantial influence on
the bifurcation locations, their variations result in negative changes to the
existence of the traveling pulse which, from a physiological perspective, may
counteract reductions in the reflection interval.

In order for pulses to propagate in excitable tissue, a balance must be struck
between the strength and timing of the local excitation and recovery
processes. The maximal strength of the currents is regulated by the
conductances and timing via the potassium channel activation and deactivation
functions. Reflections can occur when the depolarizing current is able to
overcome the repolarization level at the pulse back and initiate a retrograde
pulse. Therefore, it is no surprise that modifying the maximal calcium and
potassium conductances considerably affects both the traveling pulse and
reflection existence intervals. Furthermore, increasing the strength of
repolarization ($G_{K}$) is the most effective at making the periodic orbit
more repelling. For higher values of $G_{K}$, the strong repolarization
current inhibits reflections as a pulse must be brought very close to the
unstable periodic orbit to generate a reflection.

What is unanticipated is how altering the timing of the repolarizing potassium
current disproportionately alters the ability for reflections over traveling
waves. Parameters $u_{3a}$ and $u_{4b}$ favorably influence both the 1D spiral
stability and existence. Decreasing $u_{3a}$ shifts the potassium activation
to lower voltages, meaning the potassium channels will open and initiate the
repolarization process earlier in the action potential cycle. Increasing
$u_{4b}$ scales the deactivation function $\beta$ to slower rates. Thus, the
potassium channels remain open longer.

Moderate increases in $u_{4a}$ favorably shift $\epsilon_{*}$ to lower values
while maintaining a high $\epsilon_{SN}$. These increases correspond to
decreasing the activation rate $\alpha$ at higher voltages, but at threshold
voltages near the peak of the slow pulse (near 0 mV), the activation rates
actually increase. Therefore, the potassium channels activate more quickly
earlier in the action potential and can block perturbations propagating from
the back of the pulse. We suspect that further increases in $u_{4a}$ yields
adverse effects and decreases $\epsilon_{SN}$, and by consequence
$\epsilon_{*}$, as high values of $u_{4a}$ results in much slower activation
rates for higher voltages.

Based on our results in the Morris-Lecar model, we speculate that decreasing
the likelihood of reflections in realistic cardiac biophysical models can be
best achieved by altering the activation rates and deactivation rates of the
repolarization processes. Specifically, shifting both the activation rates and
deactivation rates to lower voltages may help most efficiently narrow the
reflection existence interval and reduce the likelihood of producing
reflections. These changes to the activation and deactivation rates act to
prolong the repolarization process by opening the gates earlier in the action
potential cycle and delaying their deactivation. Additionally, increasing the
maximum strength of the repolarization current may also have a similar
influence, but there may also be larger adverse impacts to the existence of
the necessary fast pulse.

Moreover, changes that push the system to be more type-II like may also lead
to a decrease in reflection behavior. Our results confirm that Ermentrout and
Rinzel’s hypothesis that reflections occur more robustly in systems with a
saddle-threshold in the local dynamics [21]. However, we caution that while
type II systems did tend to diminish reflection behavior, these systems also
displayed a lower range of traveling wave existence and there was not a
strictly monotonic relationship.

Furthermore, our methods described in Section 3.2 can be used to solve for any
pattern that can be formulated as time-periodic source defect and investigate
the underlying dynamics, stability, and existence. Reflected pulses are not
unique to cardiac systems; reflections and similar patterns in the wake of
traveling pulses have been observed in neurological systems [14, 27, 8, 9],
CO-dioxide reactions [39], and general pattern forming systems [41, 46].

Our work presented here is a first step in understanding the dynamical
processes that promote to the formation of reflection-induced arrhythmia and
broad physiological conditions under which they occur. Moreover, several open
questions remain. We introduced novel numerical methods to compute the
heteroclinic bifurcation point and provided numerical evidence that the
unstable periodic orbit arises through a global rearrangement of the unstable
manifold of the slow pulse, but this proposition remains to be analytically
proven. We speculate that the 1D spirals analyzed here and one-sided pulse
generators of [53] may be linked. It is likely that the two structures are
connected in an unknown way in parameter space. Finally, the Morris-Lecar is
an idealized model with links to realistic processes of cardiac dynamics, but
further investigation of the physiological processes involved in increased
reflections should be studied by applying this framework to biophysically
realistic models of cardiac tissue.

##  7 Appendix

###  7.1 Morris-Lecar Model

Ranges and default values of parameters used throughout this research are
included in Table 3. We write potassium channel dynamics of the Morris-Lecar
model using the form of activation $\alpha$ and deactivation $\beta$ rates.
These rates are related to the standard asymptotic gating variable
$n_{\infty}$ and time constant $\tau_{n}$ by

$\displaystyle n_{\infty}(V)=\frac{\alpha(V)}{\alpha(V)+\beta(V)},\ \
\tau_{n}(V)=\frac{1}{\alpha(V)+\beta(V)}.$

Table 3: Morris-Lecar system parameters, descriptions, default values, and tested parameter ranges.  Parameter | Description | Default | Range
---|---|---|---
$G_{Ca}$ | Calcium conductance | 4.4 | [4.0, 4.8]
$G_{K}$ | Potassium conductance | 8 | [6, 17]
$u_{3a}$ | Potassium activation $\alpha$ shift | 2 | [-8.5, 4.5]
$u_{3b}$ | Potassium deactivation $\beta$ shift | 2 | [-21.5, 13]
$u_{4a}$ | Potassium activation $\alpha$ scaling | 10 | [7.25, 30]
$u_{4b}$ | Potassium deactivation $\beta$ scaling | 10 | [4.5, 22]

###  7.2 Numerical Methods

We provide implementation details of the numerical methods used and developed
in this paper. All computations were performed in MATLAB 2020a [36]. Code to
compute the 1D spiral waves, heteroclinic bifurcation, and other relevant
computations can be found on GitHub [18].

####  7.2.1 Refractory-Pulse Experiment and Bisection Method

To create a functional heterogeneity, a Gaussian bump
$g(x)=B\exp\left(\frac{-(x-x_{0})^{2}}{\sigma^{2}}\right)$ is added to the
recovery variable $n(x,0)$. Time evolution is carried out with the second-
order IMEX Crank-Nicholason Adams-Bashforth splitting scheme on a periodic
domain of length 1 discretized with 1,000 fourth-order centered finite
difference spatial grid points and a temporal step size of $dt=0.03$.

The interval of amplitudes $[B_{\text{min}},B_{\text{max}}]$ that lead to at
least one reflection was computed with a series of bisection methods. First,
an amplitude that lead to reflections $B_{\text{ref}}$ was determined by
applying the bisection method to $B$ starting with values of $B=0$ and $B=1$.
Amplitudes were classified as: less than $B_{\text{min}}$ if $V(x,t_{f})$
resembled a fast pulse with one large-amplitude peak, greater than
$B_{\text{max}}$ if $V(x,t_{f})$ approached the homogeneous rest state, and in
the interval $[B_{\text{min}},B_{\text{max}}]$ if reflections were detected,
i.e. $V(x,t_{f})$ contained more than one counter-propagating large-amplitude
peak. The interval bounds $B_{\text{min}}$ and $B_{\text{max}}$ were found by
additional bisection methods between $B=0$ and $B_{\text{ref}}$, and
$B_{\text{ref}}$ and $B=1$, respectively. A stopping threshold of 1e-10 was
used in all bisection methods.

####  7.2.2 Computation of Traveling Pulses and Wave Trains

Slow and fast pulses are equilibrium solutions of

$\displaystyle 0=DU_{\xi\xi}+cU_{\xi}+F(U;\mu),\ \ \xi\in[0,L)$ (7)

in the co-moving frame $\xi=x-ct$ where $c$ is the speed of the slow or fast
pulse. Pulses are solved on a periodic discretized domain with analytical
derivatives replaced by fourth-order centered finite difference
differentiation matrices. The speed $c_{f}$ is a free variable and is solved
for with the addition of the phase condition

$\displaystyle\int_{0}^{L}\langle
U^{\prime}_{old}(\xi),U(\xi)-U_{old}(\xi)\rangle d\xi$ (8)

where the reference solution $U_{old}(\xi)$ fixes the translational symmetry
of the system and selects a unique solution [10].

Periodic wave trains are solved as stationary solutions on the co-moving
$2\pi$-periodic domain $\zeta=\kappa x-\omega t$

$\displaystyle 0=\kappa^{2}DU_{\zeta\zeta}+\omega U_{\zeta}+F(U;\mu),\ \
\zeta\in\mathcal{S}^{1}$

with wave number $\kappa$ and frequency $\omega$ connected via the nonlinear
dispersion relation $\omega=\omega(\kappa)$. A unique solution is selected by
augmenting the system with the phase condition (8). Fourier differentiation
matrices are used to approximate the derivatives.

####  7.2.3 Spectrum of Slow Pulse

To confirm the 1D spiral does not emerge through a local bifurcation of the
slow pulse, we check its essential and discrete spectrum. The discrete
spectrum is found by linearizing the reaction-diffusion equation around the
slow pulse $U_{s}(\xi)$ and solving for eigenvalues $\lambda$ of the linear
operator

$\displaystyle\mathcal{L}=D\partial_{\xi\xi}+c_{s}\partial_{\xi}+F_{U}(U_{s};\mu).$

The essential spectrum is computed as described in [42] and provides
information about the stability of the rest state at the pulse tails. As shown
in Figure 15, the essential spectra is stable for values of
$\epsilon=0.2<\epsilon_{*}$ and $\epsilon=0.35>\epsilon_{*}$. The discrete
spectrum only has a single unstable eigenvalue for these same values, which
corresponds to the one-dimensional instability generated at the saddle node
bifurcation $\epsilon_{SN}$.

Figure 15:  Essential and discrete spectrum of the slow pulse in the co-moving
frame moving at speed $c_{s}$.

####  7.2.4 Computation of 1D Spiral Waves

We consider the 1D spiral wave as a spatiotemporal equilibrium solution of

$\displaystyle 0=DU_{xx}-\omega U_{\tau}+F(U;\mu),\ \ (x,\tau)\in[-L,L]\times
S^{1}$ (9)

with $\tau=\omega t$ and angular frequency $\omega$. Following the methods of
[34, 22, 19], we split the pattern into a localized core solution $W(x,\tau)$
and the regular far-field periodic wave train pattern $U_{\infty}(\kappa
x-\tau)$ by

$\displaystyle U_{*}(x,\tau)=W(x,\tau)+\chi(x)U_{\infty}(\kappa x-\tau).$

The wave number $\kappa$ and angular frequency $\omega$ of periodic wave
trains are related via the nonlinear dispersion relation
$\omega=\omega(\kappa)$. Both $\omega$ and $\kappa$ are free variables to be
solved for. In the implementation, the far-field periodic wave train
$U_{\infty}(\kappa x-\tau)$ is first solved for on a one-dimensional
$2\pi$-periodic domain with free variable $\kappa$ using Fourier spectral
differentiation matrices. The periodic wave train is interpolated via Sinc
interpolation onto the domain $\Omega=[-L,L]\times S^{1}$ and constricted to
the far-field with the cut-off function
$\chi(x)=1-0.5\left(\tanh[m(x+d)]-\tanh[m(x-d)]\right)$ with $m=10$ and
$d=0.3$.

The core function $W(x,\tau)$ is then solved for on domain $\Omega$ with free
variable $\omega$ as the solution to (9) with Dirichlet boundary conditions
$W(x=\pm L,\tau)=0$. The phase condition

$\displaystyle
0=\int_{2\pi/\kappa}^{L}\int_{0}^{2\pi}U^{\prime}_{\infty}(\kappa
x-\tau)W(x,\tau)d\tau dx$

matches the core and far-field solutions.

Fourth-order centered finite difference differentiation matrices were used in
space, with Fourier spectral matrices in the temporal direction. The domain
size $L=1$ was used to allow for fast computations and allow the pattern to
converge to the far-field solution. On this domain size, 401 spatial grid
points and 64 temporal grid points provided more than enough accuracy for most
computations. More grid points do need to be added if the temporal period
increases rapidly (i.e., approaching the heteroclinic bifurcation) or for very
small pulse widths. An initial condition was found via numerical time
evolution and interpolated onto the spatiotemporal domain $\Omega$.

Once a 1D spiral has been computed as a spatiotemporal equilibrium solution
$U_{*}(x,\tau)$, the discrete eigenvalues can be numerically computed from the
linearized operator

$\displaystyle\mathcal{L}=D\partial_{xx}-\omega\partial_{\tau}+F_{U}(U_{*}(x,\tau);\mu)$

Leading eigenvalues were computed using the MATLAB function eigs. The domain
size and numerical discretization are the same as those for computing the
solution. The discrete eigenvalues for the default parameters with
$\epsilon=0.2$ are shown in Figure 16a. Vertically periodic branches of
eigenvalues arise due to the Floquet ambiguity in time. The form of the
leading eigenfunction, corresponding to $\lambda=2.35$, is shown in Figure
16b. The unstable eigenfunction is localized at the core region, meaning the
1D spiral pattern destabilizes at the core.

Figure 16:  (a) Eigenvalues of the 1D spiral with default parameters and
$\epsilon=0.2$. (b) Leading unstable eigenfunction, corresponding to
$\lambda=2.35$.

####  7.2.5 Computation of heteroclinic bifurcation points

Solving for the heteroclinic bifurcation point requires computing the fast
$U_{f}$ and slow $U_{s}$ pulses, the unstable eigenfunction of the slow pulse
$V_{s}^{u}$, a oppositely propagating slow-fast pulse solution $M$ on the
manifold $M_{sF}(d)$, and the unstable adjoint eigenfunction $\psi^{*}$ of
$M_{sF}(d)$. We first describe how these components are computed, then provide
the full system of equations for the heteroclinic bifurcation point.

The fast and slow pulses are solved for on a periodic one-dimensional spatial
domain as a stationary solution of (7) augmented with the standard traveling
wave phase condition (8). The unstable eigenfunction of the slow pulse
$V_{s}^{u}$ solves the eigenvalue problem linearized about the slow pulse

$\displaystyle 0$
$\displaystyle=DV_{\xi\xi}+c_{s}V_{\xi}+F_{U}(U_{s};\mu)V-\lambda V$
$\displaystyle 0$ $\displaystyle=\langle V,V\rangle-1$ (10)

where $\xi=x-c_{s}t$ is the traveling wave coordinate. Equation (10) ensures
the eigenfunction is normalized and acts as a phase condition for the
eigenvalue $\lambda\in\mathbb{C}$. The

The solution of oppositely propagating fast and slow pulses $M(\xi)$ is
constructed by combining the fast and slow pulse solutions. The peak of the
leftward traveling slow pulse is fixed at $\xi=q_{s}$, with the peak of the
rightward fast pulse set to $q_{f}=q_{s}+d$. Due to the speeds of propagation
of the pulses, we find that a spatial domain of $[-1,1]$ with $q_{s}=-0.38$
and $d=0.64$ is large enough for the fast and slow pulses to be sufficiently
separated and ensure there are no boundary effects on the computations.

The final component is the unstable adjoint eigenfunction $\psi_{*}(\xi)$,
which is solved for as a stationary solution of the adjoint eigenvalue problem

$\displaystyle 0$
$\displaystyle=D\psi_{\xi\xi}-c_{s}\psi_{\xi}+F_{U}^{*}(U_{s};\mu)\psi-\lambda\psi$
$\displaystyle 0$ $\displaystyle=\langle\psi,\psi\rangle-1$

where the second equation normalizes the adjoint eigenfunction. In computing
$\psi_{*}$, we assume that the unstable direction of the manifold $M_{sF}$
corresponds to the unstable direction of the slow pulse; thus the equation is
that of adjoint eigenvalue problem of the slow pulse.

To solve for the heteroclinic connection, the time is scaled by $T$ to
increase efficiency and we set $\tau=t/T\in[-1,1]$. The heteroclinic
connection is solved for in the co-moving frame of the slow pulse. Thus, the
full system to solve for $U(\xi,\tau)=[V(\xi,\tau),n(\xi,\tau)]^{T}$ on the
domain $[-1,1]\times[-1,1]$ is

$\displaystyle 0$ $\displaystyle=T\left[\delta
V_{\xi\xi}+c_{s}V_{\xi}+f_{1}(V,n;\mu)\right]-V_{\tau}$ (11) $\displaystyle 0$
$\displaystyle=T\left[c_{s}n_{\xi}+f_{2}(V,n;\mu)\right]-n_{\tau}$ (12)
$\displaystyle 0$
$\displaystyle=\frac{1}{T}\int_{-1}^{1}\langle\partial_{\tau}V_{old}(\bar{\xi},\tau),V(\bar{\xi},\tau)-V_{old}(\bar{\xi},\tau)\rangle
d\tau$ (13) $\displaystyle 0$
$\displaystyle=\int_{0}^{L}\langle\partial_{\xi}V_{old}(\xi,\bar{\tau}),V(\xi,\bar{\tau})-V_{old}(\xi,\bar{\tau})\rangle
d\xi$ (14) $\displaystyle 0$ $\displaystyle=\langle
V(\xi,\tau=1)-M(\xi),\psi^{*}(\xi)\rangle.$ (15)

Equations (11)-(12) solve for the $V$ and $n$ components of the heteroclinic
connection. The temporal and spatial phase conditions are given by equations
(13) - (14). The temporal phase condition acts at fixed $\xi=\bar{\xi}$ and
spatial phase condition acts at fixed $\tau=\bar{\tau}$ with reference
solution $V_{old}(\xi,\tau)$. These phase conditions are extensions of the
standard traveling wave phase condition in Equation (8). Equation (15) is the
heteroclinic boundary condition at $\tau=1$. Free parameters are $\epsilon$,
$c_{s}$, and $T$ and are additionally solved in the system and are determined
by the temporal phase condition, spatial phase condition, and heteroclinic
boundary condition, respectively.

The $\tau=-1\ (t=-T)$ boundary condition is implemented as an external
Dirichlet boundary condition. That is, the solution at $\tau=-1$ is separately
defined as $U(\tau=-1,\xi)=U_{s}(\xi)+\gamma V^{u}_{s}(\xi)$ for a small
positive $\gamma\in\mathbb{R}$, and the heteroclinic solution $U(\xi,t)$ is
solved for on the remaining temporal nodes.

The heteroclinic connection is solved for on the domain
$\Omega=[-1,1]\times[-1,1]$ using fourth-order centered-finite difference
differentiation matrices with 500 spatial grid points and 200 temporal grid
points.

## Acknowledgments

The authors express their gratitude to Björn Sandstede for invaluable
conversations on numerically computing the heteroclinic bifurcation point.

## References

  * [1] E. L. Allgower and K. Georg.  Introduction to Numerical Continuation Methods (pp 1-397).  Society for Industrial and Applied Mathematics, 2003.
  * [2] C. Antzelevitch.  Clinical applications of new concepts of parasystole, reflection, and tachycardia.  Cardiol Clin, 1:39–50, 1983.
  * [3] C. Antzelevitch and A. Burashnikov.  Overview of basic mechanisms of cardiac arrhythmia.  Card Electrophysiol Clin, 3(1):23–45, 2011.
  * [4] C. Antzelevitch and M. G. K. Jalife J.  Characteristics of reflection as a mechanism of reentrant arrhythmias and its relationship to parasystole.  Circulation, 61:182–191, 1980.
  * [5] C. Antzelevitch and G. K. Moe.  Electrotonically-mediated delayed conduction and reentry in relation to “slow responses” in mammalian ventricular conducting tissue.  Circ Res, 49:1129–39, 1981.
  * [6] U. M. Ascher, R. M. M. Mattheij, and R. D. Russel.  Numerical Solution of Boundary Value Problems for Ordinary Differential Equations.  Society for Industrial and Applied Mathematics, 1995.
  * [7] D. S. Auerbach, K. R. Grzȩda, P. B. Furspan, P. Y. Sato, S. Mironov, and J. Jalife.  Structural heterogeneity promotes triggered activity, reflection and arrhythmogenesis in cardiomyocyte monolayers.  J Physiol, 589(9):2363–2381, 2011.
  * [8] S. A. Baccus.  Synaptic facilitation by reflected action potentials: Enhancement of transmission when nerve impulses reverse direction at axon branch points.  Proceedings of the National Academy of Sciences, 95(14):8345–8350, 1998.
  * [9] S. A. Baccus, B. D. Burrell, C. L. Sahley, and K. J. Muller.  Action potential reflection and failure at axon branch points cause stepwise changes in EPSPs in a neuron essential for learning.  J Neurophysiol, 83(3):1693–1700, 2000.
  * [10] W.-J. Beyn, A. Champneys, E. Doedel, W. Govaerts, Y. A. Kuznetsov, and B. Sandstede.  Numerical Continuation and Computation of Normal Forms.  In B. Fiedler, editor, Handbook of Dynamical Systems, pages 149–219. North-Holland, 1995.
  * [11] C. Cabo and R. C. Barr.  Propagation model using the DiFrancesco-Noble equations. Comparison to reported experimental results.  Med Biol Eng Comput, 30:292–302, 1992.
  * [12] G. A. Carpenter.  Periodic Solutions of Nerve Impulse Equations.  Journal of Mathematical Analysis and Applications, 58:152–173, 1977\.
  * [13] Y. Chagnac-Amita and B. W. Connors.  Horizontal spread of synchronized activity and its control by GABA-mediated inhibition.  J Neurophys, 61:747–758, 1989.
  * [14] A. H. Cheung.  Spatiotemporal Patterns in the Wake of Traveling Waves Solutions to the Morris-Lecar Model of Neural Tissue.  PhD thesis, Boston University, 2019.
  * [15] E. N. Cytrynbaum.  Using Low Dimensional Models to Understand Cardiac Arrhythmias.  PhD thesis, University of Utah, Salt Lake City, UT, 2001.
  * [16] E. N. Cytrynbaum and T. J. Lewis.  A Global Bifurcation and the Appearance of a One-Dimensional Spiral Wave in Excitable Media.  SIAM J. Appl. Dyn. Sys., 8(1):348–370, 2009.
  * [17] J. M. Davidenko and C. Antzelevitch.  The effects of milrinone on action potential characteristics, conduction, automaticity, and reflected reentry in isolated myocardial fibers.  J Cardiovasc Pharmacol, 7:341–9, 1985.
  * [18] S. Dodson and T. J. Lewis.  GitHub Repository: 1d_spiral_waves.  https://github.com/sdodson5/1d_spiral_waves, 2021.
  * [19] S. Dodson and B. Sandstede.  Determining the Source of Period-Doubling Instabilities in Spiral Waves.  SIAM J. Appl. Dyn. Sys., 18(4):2202–2226, 2019.
  * [20] E. Doedel and L. S. Tuckerman.  Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems.  IMA Vol. Math. Appl. 119 Springer, New York, 2000.
  * [21] G. B. Ermentrout and J. Rinzel.  Reflected Waves in An Inhomogeneous Excitable Medium.  SIAM J. Appl. Math., 56(4):1107–1128, 1996.
  * [22] R. Goh and A. Scheel.  Pattern-forming fronts in a Swift-Hohenberg equation with directional quenching - parallel and oblique stripes.  J. London Math. Soc, 98(1):104–128, 2018.
  * [23] S. S. Goldstein and W. Rall.  Changes of action potential shape and velocity for changing core conductor geometry.  J. Biophys, 14:731–757, 1974.
  * [24] J. Goulet and G. B. Ermentrout.  The mechanisms for compression and reflection of cortical waves.  Biol Cybern, 105:253–268, 2011.
  * [25] S. P. Hastings.  ON THE EXISTENCE OF HOMOCLINIC AND PERIODIC ORBITS FOR THE FITZHUGH-NAGUMO EQUATIONS.  The Quarterly Journal of Mathematics, 27(1):123–134, 03 1976.
  * [26] A. J. Homburg and B. Sandstede.  Homoclinic and Heteroclinic Bifurcations in Vector Fields.  In H. Broer, F. Takens, and B. Hasselblatt, editors, Handbook of Dynamical Systems III, pages 379–524. North-Holland, 2010.
  * [27] J. F. Howe, W. H. Calvin, and J. D. Loeser.  Impulses reflected from dorsal root ganglia and from focal nerve injuries.  Brain Res., 116:139–144, 1976.
  * [28] J. Jalife and G. K. Moe.  Excitation, conduction, and reflection of impulses in isolated bovine and canine Purkinje fibers.  Circ Res, 49:233–247, 1981.
  * [29] S. M. Kandel and R. BJ.  The Mechanism of Reflection Type Reentry: A Simulation Study.  J Cardiovasc Electrophysiol, 26(12):1370–1375, 2015.
  * [30] T. Kapitula and K. Promislow.  Spectral and Dynamical Stability of Nonlinear Waves.  Springer, New York, 2013.
  * [31] N. Kopell and L. N. Howard.  Target patterns and horseshoes from a perturbed central-force problem: Some temporally periodic solutions to reaction-diffusion equations.  Studies in Applied Mathematics, 64:1–56, 1981.
  * [32] B. Krauskopf, H. M. Osinga, and J. Galan-Vioque.  Numerical Continuation Methods for Dynamical Systems.  Springer, New York, 2007.
  * [33] T. J. Lewis.  The Effects of Nonexcitable Regions on Signal Propagation in Excitable Media: Propagation Failure and Reflection.  PhD thesis, University of Utah, Salt Lake City, UT, 1998.
  * [34] D. J. B. Lloyd and A. Scheel.  Continuation and Bifurcation of Grain Boundaries in the Swift–Hohenberg Equation.  SIAM J. Appl. Dyn. Sys., 16(1):252–293, 2017.
  * [35] A. Lukas and C. Antzelevitch.  Reflected reentry, delayed conduction, and electrotonic inhibition in segmentally depressed atrial tissues.  Can J Physiol Pharmacol, 67:757–64, 1989.
  * [36] MATLAB.  version 9.8.0.1380330 (R2020a).  The MathWorks Inc., Natick, Massachusetts, 2020.
  * [37] C. Morris and H. Lecar.  Voltage oscillations in the barnacle giant muscle fiber.  Biophys J, 35(1):193–213, 1981.
  * [38] Y. Nishiura, T. Teramoto, X. Yuan, and K.-I. Ueda.  Dynamics of traveling pulses in heterogeneous media.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 17(3), 2007.
  * [39] M. Or-Guil, J. Krishnan, I. G. Kevrekidis, and M. Bär.  Pulse bifurcations and instabilities in an excitable medium: Computations in finite ring domains.  Physical Review E, 64(4):046212–9, Sept. 2001.
  * [40] A. Prat, Y.-X. Li, and P. Bressloff.  Inhomogeneity-induced bifurcation of stationary and oscillatory pulses.  Physica D, 202:177–199, 2005.
  * [41] J. D. M. Rademacher.  Homoclinic Bifurcation from Heteroclinic Cycles with Periodic Orbits and Tracefiring of Pulses.  PhD thesis, University of Minnesota, May 2004.
  * [42] J. D. M. Rademacher, B. Sandstede, and A. Scheel.  Computing absolute and essential spectra using continuation.  Physica D, 229(2):166–183, 2007.
  * [43] J. Rinzel and G. B. Ermentrout.  Analysis of neural excitability and oscillations.  In C. Koch and I. Segev, editors, Methods in Neuronal Modeling: From Synapses to Networks, pages 135–169. MIT Press, Cambridge, MA, 1989.
  * [44] G. J. Rozanski, J. Jalife, and G. K. Moe.  Reflected reentry in nonhomogeneous ventricular muscle as a mechanism of cardiac arrhythmias.  Circulation, 69:163–73, 1984.
  * [45] B. Sandstede.  Stability of travelling waves.  In B. Fiedler, editor, Handbook of Dynamical Systems II, pages 983–1055. North-Holland, 2002.
  * [46] B. Sandstede and A. Scheel.  Defects in Oscillatory Media: Toward a Classification.  SIAM J. Appl. Dyn. Sys., 3(1):1–68, 2004.
  * [47] S. H. Strogatz.  Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering. Westview Press, a member of the Perseus Books Group, 2000.
  * [48] T. Teramoto, X. Yuan, M. Bar, and Y. Nishiura.  Onset of unidirectional pulse propagation in an excitable medium with asymmetric heterogeneity.  Physical Review E, 79(4):046205, 2009.
  * [49] N. M. Van Hemel, C. A. Swenne, J. M. de Bakker, J. J. Defauw, and G. M. Guiraudon.  Epicardial reflection as a cause of incessant ventricular bigeminy.  PACE, 11:1036–1044, 1988.
  * [50] A. L. Wit, P. F. Cranefield, and B. F. Hoffman.  Slow conduction and reentry in the ventricular conducting system. ii. single and sustained circus movement in networks of canine and bovine purkinje fibers.  Circ Res, 30:11–22, 1972.
  * [51] J. D. Wright.  Separating Dissipative Pulses: The Exit Manifold.  J Dyn Diff Equ, 21(2):315–328, 2009.
  * [52] W. Xu, X. Huang, K. Takagaki, and J.-Y. Wu.  Compression and Reflection of Visually Evoked Cortical Waves.  Neuron, 55(1):119–129, 2007.
  * [53] M. Yadome, Y. Nishiura, and T. Teramoto.  Robust Pulse Generators in an Excitable Medium with Jump-Type Heterogeneity.  SIAM Journal on Applied Dynamical Systems, 13:1168–1201, 2014.
  * [54] Y. Zhou and J. Bell.  Study of propagation along nonuniform excitable fibers.  Math. Biosci., 119:169–203, 1994.