#  Annihilation dynamics during spiral defect chaos revealed by particle
models

Timothy J Tyree  Department of Physics, University of California, San Diego,
CA  Patrick Murphy  Department of Mathematics and Statistics, San Jose State
University, San Jose, CA  Wouter-Jan Rappel  Department of Physics, University
of California, San Diego, CA

###### Abstract

Pair-annihilation events are ubiquitous in a variety of spatially extended
systems and are often studied using computationally expensive simulations.
Here we develop an approach in which we simulate the pair-annihilation of
spiral wave tips in cardiac models using a computationally efficient particle
model. Spiral wave tips are represented as particles with dynamics governed by
diffusive behavior and short-ranged attraction. The parameters for diffusion
and attraction are obtained by comparing particle motion to the trajectories
of spiral wave tips in cardiac models during spiral defect chaos. The particle
model reproduces the annihilation rates of the cardiac models and can
determine the statistics of spiral wave dynamics, including its mean
termination time. We show that increasing the attraction coefficient sharply
decreases the mean termination time, making it a possible target for
pharmaceutical intervention.

††preprint: AIP/123-QED

> Many physical systems exhibit annihilation events during which pairs of
> objects collide and are removed from the system. These events occur in a
> number of soft-matter and active-matter systems that exhibit spatiotemporal
> patterning. For example, topological defects in nematic liquid crystals can
> develop motile topological defects that annihilate when they meet
> Doostmohammadi  _et al._ (2018); Liu, Shen, and Yang (2007). Pair-wise
> annihilation of defects or singularities also plays a role in a number of
> biological systems. In bacterial biofilms, for instance, imperfect cell
> alignment results in point-like defects that carry half-integer topological
> charge and can annihilate in pairs. These topological defects explain the
> formation of layers and have been proposed as a model for the buckling of
> biofilms in colonies of nematically ordered cells Boyer  _et al._ (2011);
> Copenhagen  _et al._ (2021).

##  I Introduction

In this study, we focus on the pair-wise annihilation that occurs during
spiral defect chaos in excitable systems. During spiral defect chaos, spiral
waves continuously break down to form new ones and are removed through
collisions with other spiral waves. While spiral defect chaos is present in a
variety of chemical and biological pattern-forming systems Coullet, Gil, and
Lega (1989); Rehberg, Rasenat, and Steinberg (1989); Ecke  _et al._ (1995);
Hildebrand, Bär, and Eiswirth (1995); Ouyang and Flesselles (1996); Egolf  _et
al._ (2000); Daniels and Bodenschatz (2002); Varela  _et al._ (2005); Beta
_et al._ (2006); Qu (2006); Höfer, Sherratt, and Maini (1995); Epstein and
Showalter (1996), perhaps its most studied example can be found in models of
cardiac tissue Karma (1994); Ten Tusscher  _et al._ (2004); Clayton  _et al._
(2011); Rappel (2022). These models naturally exhibit spiral waves and the
tips of these spiral waves undergo stochastic annihilation events Lilienkamp
and Parlitz (2018); Vidmar and Rappel (2019); Lilienkamp and Parlitz (2020).
Importantly, these annihilations underlie cardiac fibrillation, characterized
by unorganized electrical wave propagation in the heart Rappel (2022).
Fibrillation in the ventricles is lethal Jalife (2000) while atrial
fibrillation, the most common cardiac arrhythmia in the world with
approximately 30 million patients in 2010, is associated with increased
morbidity and mortality Naccarelli  _et al._ (2009); Miyasaka  _et al._
(2007); Chugh  _et al._ (2014); Lane  _et al._ (2017).

Previous studies have described the statistics of spiral wave tips in
spatially extended systems in the form of a stochastic birth-death process
Gil, Lega, and Meunier (1990); Daniels and Bodenschatz (2002); Beta  _et al._
(2006). More recently, this approach was applied to cardiac models, where the
creation and annihilation rates of spiral tips were determined for various
domain sizes Vidmar and Rappel (2019); Rappel  _et al._ (2022). Since $N=0$ is
an absorbing state, the dynamical state will eventually terminate. Using these
rates, the mean termination time, $\tau$, was computed and was shown to be
exponentially distributed, consistent with experiments and clinical data
Dharmaprani  _et al._ (2019); Vidmar and Rappel (2019); Rappel  _et al._
(2022). This termination time is a quantity of interest in the context of
cardiac dynamics as termination indicates the heart has transitioned into
normal sinus rhythm. Thus, minimizing $\tau$ is of critical importance for
managing cardiac fibrillation.

Previous work has shown $\tau$ depends on the tissue size, $A$, and reducing
$\tau$ can be achieved by decreasing $A$ Qu (2006); Vidmar and Rappel (2019).
Unfortunately, decreasing the effective size of cardiac tissue is not
practical and determining the dependence of $\tau$ on other physiological
properties is therefore desirable, especially if these properties can be drug-
targeted. In this study, we propose targeting the attraction coefficient, $a$,
which controls the strength of attraction between spiral wave tips.

We first simulate two spatially extended cardiac models and show that the
spiral tips that annihilate in these models display an apparent attractive
interaction at short distances and diffusive Brownian behavior at large
distances. We then formulate a stochastic particle model in which tips are
represented as particles and show that it can capture the attractive and
diffusive properties of the tips in the cardiac models. Furthermore, we show
that this particle model generates tip dynamics that reproduce both the
annihilation rates as a function of the density of tips and the distribution
of termination times for the two cardiac models. Finally, we show that
increasing the attraction coefficient significantly decreases the mean
termination time of spiral defect chaos.

##  II Cardiac Models

To determine the dynamics of spiral wave tips in the cardiac models, we
integrated the mono-domain equations, which describe the time evolution of the
membrane potential, $u$, by the excitable reaction-diffusion equation

$\partial_{t}u=D_{u}\nabla^{2}u-I_{\text{ion}}/C_{\text{m}}$ (1)

where $I_{\text{ion}}$ represents the transmembrane currents,
$C_{\text{m}}=1\mu\text{F}/\text{cm}^{2}$ is the transmembrane capacitance per
unit area, and $D_{u}$ is the diffusion coefficient Rappel (2022). To stress
the generality of our approach, we used two commonly-employed models for
cardiac tissue: the Fenton-Karma (FK) model Fenton and Karma (1998a) and the
Luo-Rudy (LR) model Luo and Rudy (1991a). The precise formulation of
$I_{\text{ion}}$ for these models is provided in Appendix A, along with model
parameters.

We integrated Eq. 1 on a square computational domain of size $A$ and enforced
periodic boundary conditions. These periodic boundary conditions resulted only
in pair-creation and pair-annihilation of spiral tips due to conservation of
vorticity in $u$ and, thus, in even numbers of tips. We used a spatial
discretization of $\Delta x=0.025$ cm and a temporal discretization of $\Delta
t=0.025$ms and used a body-centered forward-time explicit Euler method with
the Laplacian operator discretized using a five-point stencil on the square
lattice.

We chose parameter values for which both models exhibit spiral defect chaos.
In other words, a single spiral wave is unstable and will break up into
multiple spiral waves. Observations of spiral tip motion began 100 ms after
the start of the simulation at time $t=0$ so as to allow periodic boundaries
enough time to become smooth. The locations of spiral tips were determined
from the intersection points of two level sets of constant voltage ($u=0.4$
for the FK model and $u=-30$ mV for the LR model). One level set was the
equipotential line at the constant threshold voltage. The other level set was
the equipotential line at the constant threshold voltage 4 ms later in time.
Linear segments of level sets were determined by linear interpolation along
the edges that connected nodes using the marching squares method described by
Lewiner et al. Lewiner  _et al._ (2003) modified to support the periodic
boundary conditions. Intersection points were determined by solving the linear
system of equations that describes two intersecting linear segments and were
recorded as spiral tip locations. For the chosen parameter values, the width
of a planar wave was 0.41 cm for the FK model and 0.15 cm for the LR model,
computed at the same thresholds used for the detection of spiral tip
locations.

##  III Cardiac Model Results

Snapshots of typical simulations displaying a pair-wise annihilation event are
shown in Fig. 1A for the LR model (upper row) and the FK model (lower row). In
these panels, $u$ is visualized using a grayscale and clockwise and
counterclockwise tips are indicated by black and yellow symbols, respectively.
Examples of single spiral tip trajectories are shown in Fig. 1B, with the FK
trajectories plotted in blue and the LR trajectories in orange. Although
spiral tips can be both created and annihilated, we focus here on the behavior
of annihilating spiral wave tips. In our simulations, only pairs of counter-
rotating spiral waves that are connected by an activation front can
annihilate. Annihilation occurs when a depolarized region acts as a wave
block, causing the activation front to shrink in length before spontaneously
disappearing.

For both models, we examined whether tips that came within a certain distance
were able to escape and move apart or always annihilated. Specifically, we
determined the number of pairs that became separated by less than 0.1 cm and
were then able to move further apart than 0.15 cm. These pairs always
consisted of spiral waves of opposite chirality: one was rotating clockwise
while the other was rotating counterclockwise Ermakova, Pertsov, and Shnol
(1989); Kalita and Dutta (2022). Examining 355 pairs in the FK model and 99
pairs in the LR revealed that these pairs always annihilated. In the remainder
of this study, we will strictly focus on annihilation events and all
subsequent results are therefore conditioned on these dynamics.

Figure 1:  A Grayscale snapshots of $u$ showing spiral defect chaos in (top)
the LR model and (bottom) the FK model with $A=25\text{ cm}^{2}$. Indicated
are the tips of clockwise (black stars) and counterclockwise (yellow stars)
rotating spiral waves. Snapshots were taken at (left) $t\textquoteright=8$ ms,
(middle) $t\textquoteright=4$ ms, and (right) $t\textquoteright=0$ ms before
an annihilation event. Annihilation can be explained by a wave-block resulting
from a depolarized area acting as a wall to spiral tip motion. B Single tip
trajectories are shown from the FK model and the LR model (orange). The
trajectories in the LR model show fewer pivots, supporting our result that the
LR model has a larger diffusion coefficient than the FK model (see also Table
1).

To quantify the dynamics of annihilating spiral tips, we first computed the
lifetime of pairs that annihilate, $\Gamma$. For this, pair-annihilation
events were determined by ordinary nearest-neighbor particle tracking subject
to periodic boundary conditions. We computed $\Gamma$ as the time of
annihilation minus the time of creation of the younger of the pair. The
distribution of these lifetimes was approximately exponential (Fig. 2A) with
the FK model producing significantly longer-lived pairs on average (Table 1).

Figure 2:  A Probability density of lifetime $\Gamma$ of annihilating pairs in
the two cardiac models. Inset is the probability density of $\Gamma$ to
visualize the short-lived tips, revealing that the abundance of short-lived
spiral tips ($<$ 15ms) was somewhat greater for the LR model relative to the
FK model. This log-log histogram has 10 bins per decade. B Probability density
of the termination times of simulations of the two cardiac models, along with
an exponential fit (solid line).

The distribution of termination times of spiral defect chaos, i.e. the time
until all spiral tips have been removed, was also found to be fitted well with
an exponential (Fig. 2B, see also Vidmar and Rappel (2019)). This can be
understood by considering not only pair annihilation but also pair creation,
which may occur when a wave back meets a wave front. As a result of this
continuous creation and annihilation process, the number of tips, $N$, reaches
a long-lived metastable state Dykman  _et al._ (1994); Assaf and Meerson
(2010). The distribution associated with this metastable state is called the
quasi-stationary distribution and can be calculated using the master equation
for the probability $P(N,t)$ to have $N=0,2,4,8,...$ tips at time $t$ Vidmar
and Rappel (2019):

$\displaystyle\frac{dP(N,t)}{dt}$ $\displaystyle=$ $\displaystyle
W_{-2}(N+2)P(N+2,t)-W_{-2}(N)P(N,t)$
$\displaystyle\;+\;W_{+2}(N-2)P(N-2,t)-W_{+2}(N)P(N,t)\,,$

where $W_{\pm 2}$ are the annihilation and creation rates. If the number of
spiral tips reaches zero, the spiral defect chaos has terminated. In other
words, $N=0$ in the above equation is an absorbing boundary, and
$W_{+2}(0)=0$. Through repeated simulations of termination events, we
determined the distribution of termination times, which is shown for both
models in Fig. 2B.

By averaging 6,043 tip trajectories at a domain size of $A=25\text{ cm}^{2}$,
we next computed the mean squared displacement (MSD) as a function of time lag
Qian, Sheetz, and Elson (1991) (Fig. 3A). For both cardiac models, the MSD was
not significantly different from linear for long timescales with exponent
values of $1.002\pm 0.012$ (FK model for time lags from 100ms to 4000ms) and
$0.996\pm 0.019$ (LR model for time lags from 20 ms to 200 ms). Multiple other
fits using different time lag windows supported this finding (Fig. 4).

Figure 3:  A MSD of spiral tips versus temporal lag. Black lines indicate
Brownian motion. B Mean radial velocity versus inverse distance between
annihilating tips. Shaded bands represent 95% confidence intervals. Dashed
line represents $dR/dt=0$.  Figure 4:  Exponents of power laws fit to MSD
plotted as a function of lag window for the FK model (left) and the LR model
(right). The resulting exponents are visualized using a color scale and 95%
confidence intervals were determined by ordinary least squares. Gray shaded
regions indicate where the 95% confidence interval contained the slope of 1,
which corresponds to Brownian motion. The remaining regions exhibited a
statistically significant difference from Brownian behavior ($p<0.05$).
Indicated by the black dots are the fits reported in the text.

Thus, in both cardiac models, the spiral wave tips can be described as
effectively undergoing Brownian motion for long timescales. At short
timescales, however, tips did not exhibit diffusive behavior. To determine the
behavior of spiral tips at these short timescales, we analyzed the movement of
annihilating spiral waves in the simulations. This analysis revealed that
annihilation occurs when the activation front connecting the tip pair is
blocked by a nearby polarized region. The block results in a rapid shrinking
of the activation front and the removal of the pair (see Movie S1). Thus, at
short timescales, the annihilating pair of tips appear to attract, which
becomes apparent from Fig. 3B where we plot the ensemble-averaged radial
velocity $dR/dt$ as a function of the distance between the tips, $R$. This
velocity is not constant, but instead becomes more negative when $1/R$
increases (and $R$ decreases). In other words, an apparent attractive force
induced the annihilation of counter-rotating pairs of spiral waves. We should
stress that this force is only apparent and does not arise from a physical
attraction between the different spiral waves.

We next sought to quantify this attractive force. Since we were interested in
the attractive force during annihilation, we only focused on pair-wise
annihilation events in our simulations, ignoring creation events. We
identified 51,452 annihilating pairs and, for a given pair annihilating at
time $t_{f}$, we determined $R$ versus time until annihilation,
$t^{\prime}\equiv t_{f}-t\geq 0$. From this, we computed the mean squared
range (MSR) by ensemble-averaging $R^{2}$ conditioned on a given $t^{\prime}$.
The results of this ensemble-averaging are shown in Fig. 5 as solid lines,
together with the 95% confidence intervals as shaded areas. On short
timescales, the MSR demonstrates oscillations, illustrating that the apparent
attractive force has an oscillatory component. As will be discussed, this
oscillatory component is simply the manifestation of the rotation of the
spirals.

Figure 5:  MSR between annihilating tips versus time until annihilation from
simulations of the FK and LR models with shaded regions corresponding to 95%
confidence intervals. Also shown are the fits of MSR from the OPM (solid
lines) and the LPM (dashed lines).

##  IV Stochastic Particle Models

###  IV.1 Oscillatory Particle Model

To gain further insight into the annihilation dynamics of the tips using the
spatially extended cardiac models is computationally expensive, especially for
large domain sizes. Therefore, we developed a model in which the spiral wave
tips are represented by moving particles subject to an oscillatory short-range
attractive force and Brownian motion with diffusion coefficient, $D$. Such a
significant simplification of a cardiac model was earlier used to describe the
chaotic tip trajectories of a single spiral wave in the presence of
heterogeneities Lombardo and Rappel (2019).

In this oscillatory particle model (OPM), we modeled the attractive force
between two annihilating tips as inversely proportional to the distance
between them. This attractive force was assumed to consist of a constant term,
parameterized by the coefficient $a_{0}$, together with an oscillatory term,
parameterized by $a_{1}$. This oscillatory term arises from the fact that each
spiral wave rotates and that the attractive force depends on the relative
location of the two activation fronts Ermakova, Pertsov, and Shnol (1989);
Kalita and Dutta (2022). This relative location is determined by the angle of
the activation front (AAF), $\theta$, with an arbitrary, fixed axis of each
spiral wave. For example, one could choose to measure $\theta$ relative to the
positive x-axis. Measuring this AAF, however, is challenging since it requires
determining the exact location of both the tip and the activation front and is
thus prone to large discretization errors. For this reason, we will not
explicitly quantify the AAF in this study.

Based on observations of the cardiac model simulations, we will take the
frequency $\omega$ and period $T_{OPM}$ of the annihilating pair to be
identical. Furthermore, as we pointed out before, the spiral waves of an
annihilating pair have opposite chirality so that AAF of the first particle
progresses as $\theta_{1}(t)=\omega t+\theta_{1}$ and that of the second one
as $\theta_{2}(t)=-\omega t+\theta_{2}$. Then, the difference between the two
AAFs can be written as $\Delta\theta/2=\omega t+\theta_{f}$, where
$\theta_{f}=(\theta_{1}-\theta_{2})/2$. This difference between the two AAFs
is maximal and minimal once during a complete rotation of the two tips. Based
on the oscillatory nature of the MSR curves in Fig. 5, we will assume that the
attractive force depends on this difference. Taken together, the deterministic
part of the OPM for one of the particles, located at $(X_{1},Y_{1})$, is
written as

$\displaystyle\frac{dX_{1}}{dt}=\frac{X_{2}-X_{1}}{R^{2}}\Big{(}a_{0}+a_{1}\cos\big{(}\omega
t+\theta_{f}\big{)}\Big{)}$ (3)
$\displaystyle\frac{dY_{1}}{dt}=\frac{Y_{2}-Y_{1}}{R^{2}}\Big{(}a_{0}+a_{1}\cos\big{(}\omega
t+\theta_{f}\big{)}\Big{)}$ (4)

where $R=\sqrt{(X_{2}-X_{1})^{2}+(Y_{2}-Y_{1})^{2}}$ is the distance between
this particle and the other particle, located at $(X_{2},Y_{2})$. The
deterministic OPM equations for this other particle are then given by
$\frac{dX_{2}}{dt}=-\frac{dX_{1}}{dt}$,
$\frac{dY_{2}}{dt}=-\frac{dY_{1}}{dt}$.

To capture the diffusive motion of the spiral tips, we included an independent
Brownian motion term to the deterministic equations for each particle. Then,
we find that the distance between annihilating tips can be modeled by the
following Langevin equation:

$dR(t)=-\frac{2}{R(t)}\Big{(}a_{0}+a_{1}\cos\big{(}\omega
t+\theta_{f}\big{)}\Big{)}dt+\sqrt{8D}dW_{t}\,,$ (5)

where $W_{t}$ is a Wiener process at time $t$ and $D$ is estimated from the
least-square slope of MSR at large ranges ($\text{MSR}>3\text{ cm}^{2}$). Note
that the factor of two arises from the pair-wise interaction. Furthermore,
since the variance of the two two-dimensional tips add, the Wiener process is
multiplied by $\sqrt{8D}$ instead of $\sqrt{2D}$.

Using the Langevin equation, it is straightforward to derive an expression for
the MSR, $<R^{2}>$. This MSR has both a deterministic and a stochastic term:

$\displaystyle\text{MSR}_{\text{OPM}}(t^{\prime})=4\Big{(}a_{0}t^{\prime}+\frac{a_{1}}{\omega}\big{(}\sin(\omega
t^{\prime}+\theta_{f})-$
$\displaystyle\sin(\theta_{f})\big{)}+2Dt^{\prime}\Big{)}.$ (6)

This result was verified by comparing this expression to the average MSR of
10,000 statistically independent explicit simulations of Eq. 5. These
simulations can be carried out by either starting with a large value of $R$
and progressing until $R$ falls below a small threshold value or by starting
with a small value for $R$ and integrating backwards in time.

The next step was to fit Eq. IV.1 to the MSR curves obtained from the cardiac
models using simulated annealing on the last 300 ms (FK) and 100 ms (LR)
before annihilation. The fits to both cardiac models for these time intervals
are excellent (mean percent error (MPE) $<$4%), as can be seen Fig. 5. The
fitted parameters, $a_{0},\,a_{1},\,\theta_{f},\text{ and }T_{OPM}$ are
reported in Table 1, together with the aforementioned estimates for $D$. In
Table 1 we also report the period of the spiral waves in the cardiac models,
$T$, determined by computing the mean number of rotations per lifetime. A
comparison between $T_{OPM}$ and $T$ reveals that these periods match
perfectly, indicating that the oscillatory component in the MSR is due to the
rotational nature of the spiral wave.

###  IV.2 Linear Particle Model

The OPM can, in principle, be used to compute annihilation rates, which can
then be compared to results from the cardiac models. However, computing these
rates with the OPM is challenging because annihilation only occurs when the
final difference in AAF between the two particles takes a specific value.
Thus, it requires tediously tracking the AAF for each particle and enforcing
the final condition $\theta(t_{f})=\theta_{f}$ at the time of annihilation. To
circumvent this problem, we simplified the OPM to a linear particle model
(LPM) by taking the attractive force between particles as linear in $1/R$ and
dropping the oscillatory term. In this case, the attractive force is
parameterized by a single coefficient $a$ and, modeling the diffusive behavior
of spiral tips as before, the LPM reads

$dR(t)=-\frac{2a}{R}dt+\sqrt{8D}dW_{t}\,,$ (7)

which results in a MSR given by

$\text{MSR}_{\text{LPM}}(t^{\prime})=4(a+2D)t^{\prime}\,.$ (8)

To relate $a$ to the parameters of the OPM, we demanded that the MSR averaged
over the exponentially distributed lifetimes, $\langle MSR\rangle$, be equal
for both particle models. For the LPM, we can derive

$\langle
MSR_{\text{LPM}}\rangle=\int_{0}^{\infty}\Big{(}\frac{dt^{\prime}}{\Gamma}e^{-t^{\prime}/\Gamma}\Big{)}4(a+2D)t^{\prime}=4(a+2D)\Gamma\,.$
(9)

while a similar calculation for the OPM gives

$\langle
MSR_{\text{OPM}}\rangle=4\Big{(}a_{0}+a_{1}\frac{\cos(\theta_{f})-\omega\Gamma\sin(\theta_{f})}{1+(\omega\Gamma)^{2}}+2D\Big{)}\Gamma\,.$
(10)

Setting $\langle MSR_{\text{LPM}}\rangle=\langle MSR_{\text{OPM}}\rangle$
results in an analytical expression for $a$ in terms of the OPM parameters
that is independent of $D$:

$a=a_{0}+a_{1}\frac{\cos(\theta_{f})-\omega\Gamma\sin(\theta_{f})}{1+(\omega\Gamma)^{2}}\,.$
(11)

The estimates of $a$ evaluated from Eq. 11 are listed in Table 1, and
corresponding $\text{MSR}_{\text{LPM}}$ plots are shown as dashed lines in
Fig. 5. Repeating the analysis for different domain sizes, $A$, revealed that
the estimate of $a$ was largely independent of $A$ for both of the cardiac
models. An example of the MSR obtained using $A=39.0625\text{ cm}^{2}$ is
shown in Fig. 6A. For this, and all other domain sizes, we were able to fit
the MSR using the oscillatory particle model (OPM), as demonstrated by the
solid lines in Fig. 6A, which allowed us to compute the value of $a$. The
resulting value of $a$ decreased slightly as the domain size increased (Fig.
6B), suggesting that the LPM can be used to simulate spiral tip annihilation
at different domain sizes using a fixed set of model parameters determined for
a single value of $A$. Furthermore, the sum $a+2D$ also changed little as the
domain size increased (Fig. 6C). The product of this sum and the mean lifetime
is proportional to the mean squared distance (see Eq. 9). These results
indicate that the average distance between annihilating particles is
determined by local properties, which are largely insensitive to the domain
size.

Figure 6:  A MSR between annihilating tips versus time until annihilation from simulations of the FK and LR models, using a larger domain size of $A=39.0625\text{ cm}^{2}$, with shaded regions corresponding to 95% confidence intervals. The solid lines correspond to fits from the OPM while the dashed lines correspond to fits from the LPM. B Computed attraction coefficient versus domain size. C Sum of attractive and diffusive forces versus domain size. Error bars indicate 95% confidence. Table 1:  Particle properties of spiral tips from the cardiac models including parameter values corresponding to OPM and LPM.  Symbol | Fenton-Karma | Luo-Rudy
---|---|---
$\Gamma$ (ms) | 105.3 $\pm$ 1.7 | 33.4 $\pm$ 0.7
$D$ (cm2/s) | 0.115 $\pm$ 0.008 | 0.42 $\pm$ 0.14
$a_{0}$ (cm2/s) | 1.407 $\pm$ 0.016 | 4.2 $\pm$ 0.3
$a_{1}$ (cm2/s) | 1.2822 $\pm$ 0.0005 | 12.180 $\pm$ 0.012
$\theta_{f}$ (radians) | -0.541 $\pm$ 0.004 | -1.165 $\pm$ 0.003
$T_{OPM}$ (ms) | 115.94 $\pm$ 0.03 | 97.36 $\pm$ 0.12
$T$ (ms) | 115.9 $\pm$ 1.9 | 97.4 $\pm$ 0.8
$a$ (cm2/s) | 1.552 $\pm$ 0.016 | 9.3 $\pm$0.3
$r$ (cm) | 0.457$\pm$0.009 | 0.314$\pm$0.003
$\kappa$ (Hz) | 15 | 75

To compute annihilation rates in using particle dynamics instead of cardiac
model simulations, we implemented the LPM, using the obtained values of $a$
and $D$, by simulating $N$ particles with locations
$(X_{1}(t),Y_{1}(t),\dots,X_{N}(t),Y_{N}(t))$ on a square domain of size
$L^{2}$ with periodic boundary conditions. Each particle $i$ is time stepped
with time step $\Delta t=10^{-2}$ms and the $X$ coordinate of particle $i$
obeys

$X_{i}(t+\Delta
t)=X_{i}(t)+\sum_{j=1}^{N}\bigg{[}a\frac{\big{(}X_{j}(t)-X_{i}(t)\big{)}}{R_{ij}}\bigg{]}\Delta
t+\sqrt{2D\Delta t}Z\,,$

where $D>0$ is the diffusion coefficient, $Z$ is a random value with zero mean
and unit variance, and $R_{ij}$ is the distance between particle $i$ and $j$.
A similar equation can be written for $Y_{i}$. For initial conditions, we
considered $N$ uniformly distributed particles at two domain sizes
($A=25\text{ cm}^{2}$ and $A=100\text{ cm}^{2}$) and updated particle
positions every $\Delta t=0.01\text{ ms}$. Pairs of particles were removed
from the simulations, and thus annihilated, at rate $\kappa$ whenever they
were closer than a reaction range, $r$. As an estimate for this reaction
range, we chose the 25th percentile of the distribution of closest distances
between non-annihilating tips. We have verified, however, that other choices
of $r$ give similar qualitative results.

##  V Linear Particle Model Results

###  V.1 Annihilation Rates

We used the LPM to compute annihilation rates and adjusted the only free
parameter ($\kappa$) to match the annihilation rates found in the cardiac
models. This is facilitated by the fact that in our previous work we showed
that the latter, computed for different values of $A$, collapsed onto a single
curve when rescaled by $A$ Vidmar and Rappel (2019). Our results show that the
annihilation rate can be approximated by a power law,
$w_{-}=M_{-}n^{\nu_{-}}$, where $w_{-}\equiv W_{-}(N)/A$ is the rescaled
annihilation rate and $n=N/A$ is the tip density (symbols Fig. 7). The fitted
LPM annihilation rates are shown as dashed lines in Fig. 7 and resulting
parameters are listed in Table 1. These fits demonstrate that the LPM can
accurately replicate the annihilation rates of the cardiac models (MPE $<$4%).
Importantly, as was the case for $a$, these fits use the same parameter values
for both system sizes. Also note that simulations of the particle model are
much more efficient than the cardiac models, especially for large domain
sizes. Specifically, the particle model simulations use $\mathcal{O}((L/\Delta
x)^{2}-N_{avg}^{2})$ fewer operations per time step, where $N_{avg}$ is the
average number of particles. For example, for $A=100$ cm2, the speed-up
exceeded $10^{4}$-fold.

Figure 7:  Mean annihilation rate versus number density for spiral tips from
the cardiac models (symbols) and their linear particle model fits (dashed
lines).

###  V.2 Varying Attraction Strength

Although other parameters of the LPM can be varied as well, we focused on the
sensitivity of the annihilation rate to the model parameter $a$, which
controls the strength of attraction between particles. By carrying out
additional simulations of the LPM, we determined how the annihilation rate
depends on $a$ while holding the remaining parameters $D$, $r$, and $\kappa$
fixed. This rate was found to increase with increasing values of $a$, which
can be understood by realizing that larger attractive forces result in distant
particles coming closer together faster. Importantly, however, we found that
the rate was always fitted well by the power law $w_{-}=M_{-}n^{\nu_{-}}$
(Fig. 8A). For both models, we found that for increasing values of $a$ the
exponent $\nu_{-}$ became smaller (Fig. 8B) while $M_{-}$ increased (Fig. 8C).

Figure 8:  A Mean annihilation rate versus number density obtained from the
LPM using parameters corresponding to the FK model for different values of $a$
(indicated by the inset color bar). Black lines are guides to the eyes,
corresponding to power laws with exponent 4/3 (upper curve) and 2 (lower
curve). B Power law exponent as a function of $a$ computed using the LPM with
parameters corresponding to both the FK and LR model. C Corresponding power
law magnitude versus $a$. Black circles in B&C represent values of $a$
corresponding to the cardiac models. Fits considered ordinary least squares
over the interval $n\in[0.2,1]$ cm-2.

###  V.3 Statistical Properties

Using the fact that the annihilation rate in the LPM can be fitted by a power
law, we were able to compute statistical properties of spiral defect chaos in
the spatially extended cardiac models. For this, we used the creation rate in
the cardiac models, which were computed by counting the rate of creation
events conditioned on number density. As with the annihilation rates, we found
the creation rates to be captured by a power law fit, as well:
$w_{+}=M_{+}n^{\nu_{+}}$ (Fig. 9A). In the remainder, we will keep this
creation rate fixed while varying $a$, and thus the annihilation rate,
according to Figs. 8B&C. Using the power laws for creation and annihilation
rates, we constructed a transition matrix from which we computed the
distributions of termination times Newman (2010); Vidmar and Rappel (2019).
These distributions were found to be exponentially distributed for all values
of $a$ (Fig. 9B), consistent with experimental data Dharmaprani  _et al._
(2019). Furthermore, the distribution tends to smaller termination time values
for increasing $a$.

Figure 9:  A Mean creation rate (triangles) and mean annihilation rate (dots) versus number density for spiral tips from the cardiac models using different domain sizes. Dashed lines correspond to power law fits (Table 2). Black dots correspond to the mean particle density. B Probability density of termination times of the LPM for increasing values of $a$ equally spaced from $a=1\text{ cm}^{2}/\text{s}$ to $a=5\text{ cm}^{2}/\text{s}$. Parameter values correspond to the FK model and $A=25\text{ cm}^{2}$.  Table 2: Power law fits for the annihilation rates and the creation rates. Also included is the mean termination time and average number of particles.  Symbol | Fenton-Karma | Luo-Rudy
---|---|---
| Cardiac | LPM | LPM | Cardiac | LPM | LPM
$A$ (cm2) | 25 | 25 | 100 | 25 | 25 | 100
$\nu_{-}$ | 1.88$\pm$0.03 | 1.871$\pm$0.012 | 1.835$\pm$0.015 | 1.638$\pm$0.017 | 1.614$\pm$0.012 | 1.611$\pm$0.017
$M_{-}$ (Hz/cm2) | 5.6$\pm$0.3 | 5.53$\pm$0.16 | 4.75$\pm$0.17 | 16.7$\pm$0.8 | 16.9$\pm$0.7 | 12.6$\pm$0.6
$\nu_{+}$ | 0.230$\pm$0.010 | 0.230 | 0.230 | 0.715$\pm$0.010 | 0.715 | 0.715
$M_{+}$ (Hz/cm2) | 0.864$\pm$0.002 | 0.864 | 0.864 | 3.28$\pm$0.10 | 3.28 | 3.28
$\tau$ (s) | 27.8$\pm$6.5 | 25.9 | 1.68$\times 10^{9}$ | 0.74$\pm$0.06 | 0.51 | 81.6
$N_{\text{avg}}$ | 8.1$\pm$0.7 | 8.1 | 17.0 | 5.0$\pm$2.4 | 4.3 | 32.3
Figure 10:  A Average tip number as a function of $A$ computed using the
cardiac models (symbols), along with the linear prediction of Eq. S12. B
Average tip number as a function of $a$ computed using the LPM with parameter
values corresponding to $A=25\text{ cm}^{2}$. The darkened symbols correspond
to the value of $a$ representing the cardiac models. C Mean termination time
versus $A$ computed using Eq. S15 (dashed lines) and separately obtained from
the cardiac models (symbols). D Mean termination time as a function of $a$
(using parameter values for $A=25\text{ cm}^{2}$). Black circles correspond to
$a$ obtained by fitting the cardiac models.

It is now also possible to obtain approximate expressions for quantities
relevant to cardiac models and fibrillation in clinical settings. For example,
the average number of particles, $N_{\text{avg}}$ can be estimated by solving
$w_{+}(n)=w_{-}(n)$, resulting in

$\displaystyle
N_{\text{avg}}=An_{\text{avg}}=A\left(\frac{M_{+}}{M_{-}}\right)^{\frac{1}{\nu_{-}-\nu_{+}}},$
(12)

Thus, the mean particle number is predicted to increase linearly with the
domain size. Fig. 10A shows the results from Eq. 12, using parameters
corresponding to $A=25$ cm2 for all domain sizes, which is possible since $a$
is insensitive to $A$ (Fig. 6B). The estimated value using Eq. 12 agrees well
with the average number of tips from the cardiac models (symbols), especially
for larger values of $A$. Furthermore, the average number of tips computed
using Eq. 12 decreases as a function of $a$ (Fig. 10B), which can be explained
by the fact that $\nu_{-}$ decreases for increasing values of $a$ (Fig. 8B).

In addition, the mean termination time can be computed using an analytic
solution Gardiner (2004). Specifically, for paired birth-death processes, the
mean termination time $\tau$ is a function of the initial (even) number of
spiral tips, $N_{0}$:

$\displaystyle\tau(N_{0})=\sum_{k=1}^{N_{0}/2}\varphi(2k-2)\sum_{j=k}^{\infty}\frac{1}{\varphi(2j)w_{+}(2j/A)A}$
(13)

where $\varphi(2k)=\prod_{i=1}^{k}{w}_{-}(2i/A)/{w}_{+}(2i/A)$ and
$\varphi(0)\equiv 1$. Since the termination time is dominated by $\tau(2)$, we
can use the power laws for the annihilation and creation rates to derive an
approximate closed expression. For large values of $A$, the termination time
can be shown to be proportional to Vidmar and Rappel (2019)

$\tau\sim{\rm exp}\left[\frac{A}{2}\int_{2/A}^{n_{avg}}{\rm
ln}\frac{w_{+}(s)}{w_{-}(s)}ds\right]$ (14)

which becomes, after substituting the expressions for $w_{\pm}$,

$\tau\sim\left(\frac{2}{A}\right)^{\nu_{-}-\nu_{+}}e^{An_{avg}(\nu_{-}-\nu_{+})/2}\,,$
(15)

Thus, consistent with earlier work using the cardiac models Qu (2006); Vidmar
and Rappel (2019), the termination time increases exponentially with $A$.
Importantly, aside from a prefactor that is independent of $A$, this
expression provides an explicit estimate for $\tau$ for any domain size,
including ones that would be prohibitively expensive to simulate directly. In
Fig. 10C, we plot the values for $\tau$ obtained from the cardiac models as
symbols together with a fit using Eq. 14. The former ones were, of course,
only obtainable for small domain sizes while the fit can be extended to
arbitrarily large domain sizes.

Finally, using the fitted power laws (Fig. 9B), we computed $\tau$ as a
function of $a$ with the other parameters fixed to those corresponding to the
cardiac models. We found that $\tau$ decreased by a factor of $\sim 10^{2}$ in
response to $a$ increasing by a mere factor of $10$, as shown in Fig. 10D.
This great sensitivity of $\tau$ to $a$, considered together with the relative
insensitivity of $a$ to $A$, suggests that changing $a$ is a mechanism for
controlling the mean termination time for minimally varying domain sizes.
Specifically, if modifying the electrophysiological parameters of the cardiac
models increases attraction only, then $\tau$ decreases, resulting in shorter
termination times.

To verify that it is possible to vary $a$, and thus $\tau$, we carried out
additional simulations of the cardiac models. For this, we altered the
excitability parameter $\tau_{d}$ in the FK model and the extracellular
potassium concentration $[K^{+}]_{o}$ in the LR model. We followed the same
procedure as before and fitted the MSR expression for the OPM to the MSR
curves obtained from the cardiac models using simulated annealing on the last
300 ms (FK) and 100 ms (LR) before annihilation. The fit results are shown in
Figure 11 for both models. These fits, together with measured values of $D$
and $\Gamma$, were then used to compute $a$ and the LPM model was simulated to
determine $\tau$.

For the FK model, we increased $\tau_{d}$ by 20% and obtained as fitted
parameter values $a_{0}=1.6077$ cm${}^{2}/$s, $a_{1}=1.3311$ cm${}^{2}/$s,
$T_{OPM}=110.3665$ ms, and $\theta_{f}=-1.0730$ radians (RMSE=0.0459 cm2). The
measured parameters were $D=0.048\pm 0.007$ cm${}^{2}/$s and
$\Gamma=101.3\pm$1.3 ms, resulting in $a=1.82\pm 0.01$ cm${}^{2}/$s, which was
17.5% larger than the corresponding value for $a$ in Table 1. This led to a
decrease of the mean termination time by 35%. In the LR model, we reduced
$[K^{+}]_{o}$ by 29.6% and found the following fitted parameter values:
$a_{0}=4.1932$ cm${}^{2}/$s, $a_{1}=5.9593$ cm${}^{2}/$s, $T_{OPM}=82.4087$
ms, and $\theta_{f}=-0.9490$ radians (RMSE=0.0352 cm2). Measured parameters
were $D=0.42\pm 0.03$ cm${}^{2}/$s and $\Gamma=45.9\pm 1.1$ ms, resulting in
$a=5.73\pm 0.07$ cm${}^{2}/$s, which was $38.6\%$ smaller than the
corresponding value for $a$ in Table 1. The termination time in this case
increased by $173.6\%$. Thus, increasing $a$ reduced $\tau$ for the FK model,
and reducing $a$ increased $\tau$ for the LR model.

Figure 11:  A MSR from the FK model with the excitability parameter,
$\tau_{d}$, increased by 20%. The blue shaded region represents 95% confidence
intervals estimated via bootstrap. The black line represents the fit to the
OPM. B MSR from the LR model with the extracellular potassium concentration,
$[K^{+}]_{o}$, decreased by 29.6%. The orange shaded region represents 95%
confidence intervals estimated via bootstrap. The red line represents the fit
to the OPM.

##  VI Concluding Remarks

In summary, this study reveals that the annihilation dynamics of spiral waves
in spatially extended cardiac models can be captured by a computationally
efficient particle model. In this model, the spiral wave tips are represented
by diffusing particles subject to a locally attractive force. We showed that
the particle model accurately reproduced the annihilation rates of spiral tips
in the cardiac models and that it can be used to efficiently model their tip
dynamics. We also used the particle model to show that increasing the strength
of the apparent attraction force accelerates annihilation, thus decreasing the
mean termination time.

What have we ignored and how can this study be extended? First of all, we
exclusively focused on annihilation dynamics. Clearly, to construct more
realistic particle models, we will have to include creation events. These
events may require the inclusion of repulsion of spiral wave tips, as has been
reported in some cardiac models Ermakova, Pertsov, and Shnol (1989). Once
these events are incorporated into the model, we should be able to determine
how parameters that describe these creation events can change termination
times. Second, we have simulated the dynamics on a domain with periodic
boundary conditions, guaranteeing that annihilation is a pair-wise event.
Introducing non-conducting boundaries will necessitate a description of the
annihilation of a single spiral wave that collides with one of these
boundaries. Third, we have simplified the oscillatory model to a linear model
such that we did not have to compute the AAFs for the spiral waves. It would
be interesting to extend our study to explicitly enforce conditions on the
relative AAFs at the time of annihilation. Finally, future work involving
cardiac models can include a more thorough investigation of the dependence of
the attraction coefficient between spiral tips on trans-membrane currents.
Reproducing desirable effects on these currents can then be the target of drug
discovery, potentially opening a new door to noninvasive therapies for
clinically significant symptoms of cardiac fibrillation.

##  VII Acknowledgements

This work was supported by the National Institute of Health through awards
HL083359, HL122384, HL149134, and T32GM127235. T.T. acknowledges the Open
Science Grid Pordes  _et al._ (2007); Sfiligoi  _et al._ (2009) for allocation
of computational resources. P.M acknowledges the support of the Research
Training in Mathematical and Computational Biology Grant (1148230) and the NSF
Division of Mathematical Sciences award 1903275.

##  Appendix A Detailed Equations and Parameters of the Full Cardiac Models

###  A.1 Fenton-Karma Model

The ionic currents in the Fenton-Karma (FK) cardiac model Fenton and Karma
(1998b) are written as:

$I_{\text{ion}}=I_{\text{fi}}+I_{\text{so}}+I_{\text{si}},$

where the fast inward current is

$I_{\text{fi}}(u,v)=-\frac{v}{\tau_{d}}(1-u)(u-u_{c})\Theta(u-u_{c}),$

the slow outward current is

$I_{\text{so}}(u)=u\frac{\Theta(u_{c}-u)}{\tau_{0}}+\frac{\Theta(u-u_{c})}{\tau_{r}},$

and the slow inward current is

$I_{\text{si}}(u,w)=\frac{w}{2\tau_{\text{si}}}\Big{(}1+\tanh\big{(}k(u-u_{c}^{\text{si}})\big{)}\Big{)}.$

In these equations, $v$ and $w$ represent the fast and slow variable,
respectively:

$\frac{dv}{dt}=(1-v)\Theta(u_{c}-u)\bigg{(}\frac{\Theta(u-u_{v})}{\tau_{v1}^{-}}+\frac{\Theta(u_{v}-u)}{\tau_{v2}^{-}}\bigg{)}-v\frac{\Theta(u-u_{c})}{\tau_{v}^{+}}$
(16)

and

$\frac{dw}{dt}=(1-w)\frac{\Theta(u_{c}-u)}{\tau_{w}^{-}}-w\frac{\Theta(u-u_{c})}{\tau_{w}^{+}}.$
(17)

Parameter values for the model were chosen to be $\tau_{d}=0.45$ms,
$\tau_{0}=12.5$ms, $\tau_{r}=$33.25ms, $\tau_{\text{si}}=29$ms, $k=10$,
$u_{c}^{\text{si}}=0.85$, $u_{c}=0.13$, $u_{v}=0.04$, $\tau_{v1}^{-}=1250$ms,
$\tau_{v2}^{-}=19.6$ms, $\tau_{v}^{+}=13.03$ms, $\tau_{w}^{-}=40$ms, and
$\tau_{w}^{+}=800$ms. Using a diffusion coefficient of $D_{u}=0.0005$cm2/ms
resulted in a conduction velocity of a planar wave of $c_{v}=51$cm/s, which is
within the electrophysiological range.

###  A.2 Luo-Rudy Model

The total current for the Luo-Rudy (LR) cardiac model Luo and Rudy (1991b)
used in this study is given by

$I_{\text{ion}}=I_{\text{Na}}+I_{\text{si}}+I_{\text{K1T}}+I_{\text{K}},$

where the first term is the sodium current

$I_{\text{Na}}=(u-E_{\text{Na}})G_{\text{Na}}jhm^{3}.$

with parameter values $G_{\text{Na}}=16\text{mS}/\text{cm}^{2}$ and
$E_{\text{Na}}=54.4\text{mV}$. The second term is the slow inward current

$I_{\text{si}}=G_{\text{si}}(u-E_{\text{si}})fd,$

where $G_{\text{si}}=0.052\text{mS}/\text{cm}^{2}$ and
$E_{\text{si}}=-82.3\text{mV}-13.0287\log\big{(}[\text{Ca}^{+2}]_{i})\big{)}$.
The last two terms represent the potassium current, which is computed using
the equations and parameters of Qu et al. in Ref. Qu, Weiss, and Garfinkel
(2000); Qu  _et al._ (2000). In these expression, the intracellular calcium is
given by

$d[\text{Ca}^{+2}]_{i}/dt=-10^{-7}I_{\text{si}}+\big{(}0.07\big{)}\big{(}10^{-7}-[\text{Ca}^{+2}]_{i}\big{)}$

Furthermore, the dimensionless gating variables, $y=m,h,j,d,f$, are described
by equations of the form

$dy/dt=(y_{\infty}-y)/\tau_{y}.$ (18)

where

$y_{\infty}(u)=a_{y}(u)\tau_{y}(u)$

and

$\tau_{y}(u)=\frac{1}{a_{y}(u)+b_{y}(u)},$

In these equations, $a_{y}$, $b_{y}$ are dimensionless monotonic functions,
which were evaluated in constant time using a lookup table based on the
equations of Qu et al. Qu, Weiss, and Garfinkel (2000); Qu  _et al._ (2000).
Additional parameters can be found in these references. For complete
reproducibility, the equations and parameters that we used are included in
Table 3 and here. We have

$\displaystyle
I_{\text{K1T}}=I_{K1}(V_{\text{m}})+(1\mu\text{A}/\text{cm}^{2})\frac{(V_{\text{m}}+87.95\text{mV})/54.6448\text{mV}}{1+e^{(7.488\text{mV}\,-V_{\text{m}})/5.98\text{mV}}}+(1\mu\text{A}/\text{cm}^{2})\frac{V_{\text{m}}+59.87\text{mV}}{25.5037\text{mV}},$

where

$\displaystyle
I_{K1}(V_{\text{m}})=\frac{(1\mu\text{A}/\text{cm}^{2})\frac{V_{\text{m}}+87.95\text{mV}}{1.62129\text{mV}}{\sqrt{\frac{[K^{+}]_{o}}{5.4\text{mM}}}}}{\left(1+e^{\frac{V_{\text{m}}+28.735\text{mV}}{4.19287\text{mV}}}\right)\left(\frac{e^{\frac{V_{\text{m}}-506.36\text{mV}}{16.1943\text{mV}}}+0.49124e^{\frac{V_{\text{m}}+93.426\text{mV}}{12.4502\text{mV}}}}{1+e^{-\frac{V_{\text{m}}+92.703\text{mV}}{1.94439\text{mV}}}}+\frac{1.02\text{mV}}{1+e^{\frac{V_{\text{m}}+28.735\text{mV}}{4.19287\text{mV}}}}\right)},$

and where we have taken $[K^{+}]_{o}=5.4\text{mM}$. We computed $I_{\text{K}}$
using $I_{\text{K}}=I_{1}x$, where

$\displaystyle
I_{1}=\Big{(}2.837G_{\text{K}}\sqrt{\frac{5.4\text{mM}}{[K^{+}]_{o}}}\Big{)}\frac{V_{\text{m}}-E_{1}}{V_{\text{m}}+77\text{mV}}\Big{(}\exp\Big{(}\frac{V_{\text{m}}+77\text{mV}}{25\text{mV}}\Big{)}-1\Big{)}\Big{/}\exp\Big{(}\frac{V_{\text{m}}+35\text{mV}}{25\text{mV}}\Big{)}\Big{)}$

where have used the values $G_{\text{K}}=0.423\text{mS}/\text{cm}^{2}$ and

$E_{1}=(10^{3}\text{mV})\frac{RT}{F}\log\Big{(}\frac{[K^{+}]_{o}+0.01833[Na^{+}]_{o}}{[K^{+}]_{i}+0.01833[Na^{+}]_{i}}\Big{)}\approx-77.61\text{mV}.$

where we took $[K^{+}]_{i}=145\text{mM}$, $[{Na^{+}]_{i}=18\text{mM}}$,
$[Na^{+}]_{o}=145\text{mM}$, $R=8.3145\text{J}/(\text{mol}\,^{\circ}\text{K})$
as the universal gas constant, and $F=96485.3321233100184\text{C}/\text{mol}$
as Faraday’s constant. We supposed a homeostatic body temperature of
$T=37^{\circ}C$ fixed. We have similarly approximated

$E_{\text{K}1}=(10^{3}\text{mV})\frac{RT}{F}\log\bigg{(}\frac{[K^{+}]_{o}}{[K^{+}]_{i}}\bigg{)}\approx
87.95\text{mV}.$

To avoid numerical overflow when time evolving these gating variables, we
treated time scales as zero, $\tau_{y}=0$, whenever they took a value
$\tau_{y}\leq 5\cdot 10^{-4}$ms. Additionally, we evaluated $I_{1}$ using
L’Hospital’s Rule in the $10^{-6}\text{mV}$ neighborhood of
$V_{\text{m}}=-77\text{mV}$.

Finally, for the diffusion coefficient we used $D_{u}=0.001$cm2/ms, resulting
in a conduction velocity of $c_{v}=33$cm/s, which is within the
electrophysiological range.

Table 3: Parameters of the Luo-Rudy-I model used in our study. For ease of notation, we used the dimensionless transmembrane voltage field $u=V_{m}/(1mV)$. Gating variables time evolve according to Eq. S18. $y$ | $a_{y}$ | $b_{y}$
---|---|---
m | $\frac{u+47.13}{3.125}\Big{/}\Big{(}1-e^{-\frac{u+47.13}{10}}\Big{)}$ | $0.08e^{-u/11}$
h | $0.135\exp\Big{(}\frac{u+80}{6.8}\Big{)}\Theta(40-u)$ |  $\frac{100}{13}\frac{\Theta(40-u)}{1+\exp\big{(}\frac{10.66-u}{11.1}\big{)}}\;$ $+\;\Big{(}3.56e^{0.079u}+3.1\cdot 10^{5}e^{0.35u}\Big{)}\Theta(u-40)$
j | $-(u+37.78)\frac{\left(3.474\cdot 10^{-5}e^{-0.04391u}+127140e^{0.24444u}\right)}{1+e^{0.311(u+79.23)}}\Theta(40-u)$ |  $\frac{0.3e^{-2.535\cdot 10^{-7}u}}{1+e^{-0.1(u+32)}}\Theta(u-40)$ \+ $0.1212e^{-0.01052u}\frac{\Theta(40-u)}{1+e^{-0.1378(u+40.14)}}$
d | $0.095e^{0.01(5-u)}\Big{/}\Big{(}1+e^{0.072(5-u)}\Big{)}$ | $0.07e^{-0.017(44+u)}\Big{/}\Big{(}1+e^{0.05(44+u)}\Big{)}$
f | $0.012e^{-0.008(28+u)}\Big{/}\Big{(}1+e^{0.15(28+u)}\Big{)}$ | $0.0065e^{-0.02(30+u)}\Big{/}\Big{(}1.\,+e^{-0.2(30+u)}\Big{)}$
x | $0.5\cdot 10^{-3}e^{0.083(u+50)}\Big{/}\Big{(}1+e^{0.057(u+50)}\Big{)}$ | $1.3\cdot 10^{-3}\exp\big{(}\frac{u+20}{25}\big{)}\Big{/}\Big{(}1+\exp\big{(}\frac{3}{2}\frac{u+20}{25}\big{)}\Big{)}$

## References

  * Doostmohammadi  _et al._ (2018) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, “Active nematics,” Nature communications 9, 1–13 (2018).
  * Liu, Shen, and Yang (2007) C. Liu, J. Shen, and X. Yang, “Dynamics of defect motion in nematic liquid crystal flow: modeling and numerical simulation,” Commun. Comput. Phys 2, 1184–1198 (2007).
  * Boyer  _et al._ (2011) D. Boyer, W. Mather, O. Mondragón-Palomino, S. Orozco-Fuentes, T. Danino, J. Hasty, and L. S. Tsimring, “Buckling instability in ordered bacterial colonies,” Physical biology 8, 026008 (2011).
  * Copenhagen  _et al._ (2021) K. Copenhagen, R. Alert, N. S. Wingreen, and J. W. Shaevitz, “Topological defects promote layer formation in myxococcus xanthus colonies,” Nature Physics 17, 211–215 (2021).
  * Coullet, Gil, and Lega (1989) P. Coullet, L. Gil, and J. Lega, “Defect-mediated turbulence,” Physical Review Letters 62, 1619 (1989).
  * Rehberg, Rasenat, and Steinberg (1989) I. Rehberg, S. Rasenat, and V. Steinberg, “Traveling waves and defect-initiated turbulence in electroconvecting nematics,” Physical Review Letters 62, 756 (1989).
  * Ecke  _et al._ (1995) R. E. Ecke, Y. Hu, R. Mainieri, and G. Ahlers, “Excitation of spirals and chiral symmetry breaking in rayleigh-bénard convection,” Science 269, 1704–1707 (1995).
  * Hildebrand, Bär, and Eiswirth (1995) M. Hildebrand, M. Bär, and M. Eiswirth, “Statistics of topological defects and spatiotemporal chaos in a reaction-diffusion system,” Physical Review Letters 75, 1503 (1995).
  * Ouyang and Flesselles (1996) Q. Ouyang and J.-M. Flesselles, “Transition from spirals to defect turbulence driven by a convective instability,” Nature 379, 143 (1996).
  * Egolf  _et al._ (2000) D. A. Egolf, I. V. Melnikov, W. Pesch, and R. E. Ecke, “Mechanisms of extensive spatiotemporal chaos in rayleigh–bénard convection,” Nature 404, 733–736 (2000).
  * Daniels and Bodenschatz (2002) K. E. Daniels and E. Bodenschatz, “Defect turbulence in inclined layer convection,” Physical Review Letters 88, 034501 (2002).
  * Varela  _et al._ (2005) H. Varela, C. Beta, A. Bonnefont, and K. Krischer, “Transitions to electrochemical turbulence,” Physical Review Letters 94, 174104 (2005).
  * Beta  _et al._ (2006) C. Beta, A. S. Mikhailov, H. H. Rotermund, and G. Ertl, “Defect-mediated turbulence in a catalytic surface reaction,” EPL (Europhysics Letters) 75, 868 (2006).
  * Qu (2006) Z. Qu, “Critical mass hypothesis revisited: role of dynamical wave stability in spontaneous termination of cardiac fibrillation,” American Journal of Physiology-Heart and Circulatory Physiology 290, H255–H263 (2006).
  * Höfer, Sherratt, and Maini (1995) T. Höfer, J. A. Sherratt, and P. K. Maini, “Dictyostelium discoideum: cellular self-organization in an excitable biological medium,” Proceedings of the Royal Society of London. Series B: Biological Sciences 259, 249–257 (1995).
  * Epstein and Showalter (1996) I. R. Epstein and K. Showalter, “Nonlinear chemical dynamics: oscillations, patterns, and chaos,” The Journal of Physical Chemistry 100, 13132–13147 (1996).
  * Karma (1994) A. Karma, “Electrical alternans and spiral wave breakup in cardiac tissue,” Chaos: An Interdisciplinary Journal of Nonlinear Science 4, 461–472 (1994).
  * Ten Tusscher  _et al._ (2004) K. Ten Tusscher, D. Noble, P. Noble, and A. Panfilov, “A model for human ventricular tissue,” American Journal of Physiology-Heart and Circulatory Physiology 286, H1573–H1589 (2004).
  * Clayton  _et al._ (2011) R. Clayton, O. Bernus, E. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V. Panfilov, F. B. Sachse, G. Seemann, and H. Zhang, “Models of cardiac tissue electrophysiology: progress, challenges and open questions,” Progress in biophysics and molecular biology 104, 22–48 (2011).
  * Rappel (2022) W.-J. Rappel, “The physics of heart rhythm disorders,” Physics Reports 978, 1–45 (2022).
  * Lilienkamp and Parlitz (2018) T. Lilienkamp and U. Parlitz, “Terminal transient phase of chaotic transients,” Physical review letters 120, 094101 (2018).
  * Vidmar and Rappel (2019) D. Vidmar and W.-J. Rappel, “Extinction dynamics of spiral defect chaos,” Physical Review E 99, 012407 (2019).
  * Lilienkamp and Parlitz (2020) T. Lilienkamp and U. Parlitz, “Terminating transient chaos in spatially extended systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 051108 (2020).
  * Jalife (2000) J. Jalife, “Ventricular fibrillation: mechanisms of initiation and maintenance,” Annual review of physiology 62, 25–50 (2000).
  * Naccarelli  _et al._ (2009) G. V. Naccarelli, H. Varker, J. Lin, and K. L. Schulman, “Increasing prevalence of atrial fibrillation and flutter in the United States,” Am. J. Cardiol. 104, 1534–1539 (2009).
  * Miyasaka  _et al._ (2007) Y. Miyasaka, M. E. Barnes, K. R. Bailey, S. S. Cha, B. J. Gersh, J. B. Seward, and T. S. Tsang, “Mortality trends in patients diagnosed with first atrial fibrillation: a 21-year community-based study,” Journal of the American College of Cardiology 49, 986–992 (2007).
  * Chugh  _et al._ (2014) S. S. Chugh, R. Havmoeller, K. Narayanan, D. Singh, M. Rienstra, E. J. Benjamin, R. F. Gillum, Y.-H. Kim, J. H. McAnulty Jr, Z.-J. Zheng, _et al._ , “Worldwide epidemiology of atrial fibrillation: a global burden of disease 2010 study,” Circulation 129, 837–847 (2014).
  * Lane  _et al._ (2017) D. A. Lane, F. Skjøth, G. Y. Lip, T. B. Larsen, and D. Kotecha, “Temporal trends in incidence, prevalence, and mortality of atrial fibrillation in primary care,” Journal of the American Heart Association 6, e005155 (2017).
  * Gil, Lega, and Meunier (1990) L. Gil, J. Lega, and J. Meunier, “Statistical properties of defect-mediated turbulence,” Physical Review A 41, 1138 (1990).
  * Rappel  _et al._ (2022) W.-J. Rappel, D. E. Krummen, T. Baykaner, J. Zaman, A. Donsky, V. Swarup, J. M. Miller, and S. M. Narayan, “Stochastic termination of spiral wave dynamics in cardiac tissue,” Frontiers in Network Physiology , 2 (2022).
  * Dharmaprani  _et al._ (2019) D. Dharmaprani, M. Schopp, P. Kuklik, D. Chapman, A. Lahiri, L. Dykes, F. Xiong, M. Aguilar, B. Strauss, L. Mitchell, _et al._ , “Renewal theory as a universal quantitative framework to characterize phase singularity regeneration in mammalian cardiac fibrillation,” Circulation: Arrhythmia and Electrophysiology 12, e007569 (2019).
  * Fenton and Karma (1998a) F. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation,” Chaos 8, 20–47 (1998a).
  * Luo and Rudy (1991a) C. H. Luo and Y. Rudy, “A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction,” Circ Res 68, 1501–1526 (1991a).
  * Lewiner  _et al._ (2003) T. Lewiner, H. Lopes, A. W. Vieira, and G. Tavares, “Efficient implementation of marching cubes’ cases with topological guarantees,” Journal of graphics tools 8, 1–15 (2003).
  * Ermakova, Pertsov, and Shnol (1989) E. Ermakova, A. Pertsov, and E. Shnol, “On the interaction of vortices in two-dimensional active media,” Physica D: Nonlinear Phenomena 40, 185–195 (1989).
  * Kalita and Dutta (2022) H. Kalita and S. Dutta, “Interaction of multiple spiral rotors in a reaction-diffusion system,” Physical Review E 105, 054213 (2022).
  * Dykman  _et al._ (1994) M. Dykman, E. Mori, J. Ross, and P. Hunt, “Large fluctuations and optimal paths in chemical kinetics,” The Journal of chemical physics 100, 5735–5750 (1994).
  * Assaf and Meerson (2010) M. Assaf and B. Meerson, “Extinction of metastable stochastic populations,” Physical Review E 81, 021116 (2010).
  * Qian, Sheetz, and Elson (1991) H. Qian, M. P. Sheetz, and E. L. Elson, “Single particle tracking. analysis of diffusion and flow in two-dimensional systems,” Biophysical journal 60, 910–921 (1991).
  * Lombardo and Rappel (2019) D. M. Lombardo and W.-J. Rappel, “Chaotic tip trajectories of a single spiral wave in the presence of heterogeneities,” Physical Review E 99, 062409 (2019).
  * Newman (2010) M. Newman,  _Networks; an introduction_ (Oxford University, Oxford, 2010).
  * Gardiner (2004) C. Gardiner,  _Handbook of stochastic methods_ , Vol. 3 (Springer Berlin, 2004) p. 260.
  * Pordes  _et al._ (2007) R. Pordes, D. Petravick, B. Kramer, D. Olson, M. Livny, A. Roy, P. Avery, K. Blackburn, T. Wenaus, F. Würthwein, I. Foster, R. Gardner, M. Wilde, A. Blatecky, J. McGee, and R. Quick, “The open science grid,” in  _J. Phys. Conf. Ser._, 78, Vol. 78 (2007) p. 012057.
  * Sfiligoi  _et al._ (2009) I. Sfiligoi, D. C. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, and F. Würthwein, “The pilot way to grid resources using glideinwms,” in  _2009 WRI World Congress on Computer Science and Information Engineering_, 2, Vol. 2 (2009) pp. 428–432.
  * Fenton and Karma (1998b) F. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation,” Chaos: An Interdisciplinary Journal of Nonlinear Science 8, 20–47 (1998b).
  * Luo and Rudy (1991b) C.-h. Luo and Y. Rudy, “A model of the ventricular cardiac action potential. depolarization, repolarization, and their interaction.” Circulation research 68, 1501–1526 (1991b).
  * Qu, Weiss, and Garfinkel (2000) Z. Qu, J. N. Weiss, and A. Garfinkel, “From local to global spatiotemporal chaos in a cardiac tissue model,” Physical Review E 61, 727 (2000).
  * Qu  _et al._ (2000) Z. Qu, A. Garfinkel, P.-S. Chen, and J. N. Weiss, “Mechanisms of discordant alternans and induction of reentry in simulated cardiac tissue,” Circulation 102, 1664–1670 (2000).