Dynamical mechanism of atrial fibrillation: a topological approach

TL;DR

This paper investigates the mechanisms sustaining atrial fibrillation using a topological approach, revealing a balance between wave coalescence and collapse, challenging the traditional view of wave breakup as the main factor.

Contribution

It introduces a novel topological framework to explain atrial fibrillation maintenance, emphasizing wave coalescence and collapse over wave breakup.

Findings

Spatiotemporal chaos maintained by wave coalescence and collapse.

Wave breakup does not sustain atrial fibrillation.

Topological processes explain the dynamics of wavelet number changes.

Abstract

While spiral wave breakup has been implicated in the emergence of atrial fibrillation, its role in maintaining this complex type of cardiac arrhythmia is less clear. We used the Karma model of cardiac excitation to investigate the dynamical mechanisms that sustain atrial fibrillation once it has been established. The results of our numerical study show that spatiotemporally chaotic dynamics in this regime can be described as a dynamical equilibrium between topologically distinct types of transitions that increase or decrease the number of wavelets, in general agreement with the multiple wavelets hypothesis. Surprisingly, we found that the process of continuous excitation waves breaking up into discontinuous pieces plays no role whatsoever in maintaining spatiotemporal complexity. Instead this complexity is maintained as a dynamical balance between wave coalescence -- a unique, previously unidentified, topological process that increases the number of wavelets -- and wave collapse -- a different topological process that decreases their number.

Equations

$${\partial }_t\mathbf{u}=D{\nabla }^2\mathbf{u}+\mathbf{f}(\mathbf{u}),$$

$$f_1$$

$$f_2$$

$$R=\{\mathbf{x}|{\partial }_t{u}_1(t,\mathbf{x}))<0\},$$

$${\dot{u}}_1$$

$${\dot{u}}_2$$

$$E=\{\mathbf{x}|f_1^{+}(\mathbf{u}(t,\mathbf{x}))=O(ϵ)\},$$

$$g(\mathbf{u})=0,$$

$$\partial E=\{\mathbf{x}|g(\mathbf{u}(t,\mathbf{x}))=0\}.$$

$$g(\mathbf{u})=\{\begin{array}{cc}{u}_1-u_1^{<},& u_2\le u_2^{<},\\ f_1^0(\mathbf{u}),& u_2^{<}<u_2<u_2^{>},\\ u_1-u_1^{>},& u_2\ge u_2^{>}.\end{array}$$

$$c=-\frac{\mathbf{b}\cdot {\partial }_t\mathbf{u}}{\mathbf{b}\cdot {\partial }_n\mathbf{u}},$$

$$E=\{\mathbf{x}|g(\mathbf{u}(t,\mathbf{x}))>0\},$$

$$g(\mathbf{u})=\mathbf{a}\cdot {\partial }_t\mathbf{u}=0,$$

$$\partial E=\{\mathbf{x}|{\partial }_t{u}_2(t,\mathbf{x})=0\}.$$

$$u_1(t,\mathbf{x})={\bar{u}}_1,f_1^0({\bar{u}}_1,u_2)=0,$$

$$u_1(t,\mathbf{x})={\bar{u}}_1,{\partial }_t{u}_1(t,\mathbf{x})=0,$$

$${\partial }_t\mathbf{u}(t,\mathbf{x})=\omega {\partial }_{\theta }\mathbf{u}(t,\mathbf{x}),$$

$${\partial }_t\mathbf{u}(t,{\mathbf{x}}_o)=0.$$

$$\partial R=\{\mathbf{x}|{\partial }_t{u}_1(t,\mathbf{x}))=0\}$$

$$\phi =\arg ({\partial }_t{u}_1+i{\partial }_t{u}_2),$$

$${\int }_{u_1^{-}}^{u_1^{+}}\mathrm{d}{u}_1{f}_1(u_1,{\bar{u}}_2)=0,$$

$$q_j=\mathrm{sign}(\hat{\mathbf{z}}\cdot \nabla u_1\times \nabla u_2).$$

$$w_j(\mathbf{x})=e^{-r_j/d_j},$$

$$J({\tilde{q}}_j)={\int }_{\Omega }{\mathrm{d}}^2\mathbf{x}{w}_j(\mathbf{x}){\left|{\partial }_t\mathbf{u}(t,\mathbf{x})-{\tilde{q}}_j\omega {\partial }_{\theta }\mathbf{u}(t,\mathbf{x})\right|}^2,$$

$$q=\sum _j{q}_j=0$$

$$D_{11}{u}_1^{\prime \prime }+c{u}_1^{\prime }+f_1(\mathbf{u})=0,$$

$${\partial }_t{u}_1=D_{11}{\nabla }^2{u}_1+f_1(\mathbf{u})=0$$

$$c=c_0-D_{11}\kappa ,$$

$$D_{11}[{\partial }_r^2{u}_1+r^{-1}{\partial }_r{u}_1+r^{-2}{\partial }_{\theta }^2{u}_1]+f_1(\mathbf{u})=0$$

$$D_{22}[{\partial }_r^2{u}_2+r^{-1}{\partial }_r{u}_2+r^{-2}{\partial }_{\theta }^2{u}_2]+f_2(\mathbf{u})=0$$

$$f_2(u_1(r_c),{\bar{u}}_2)=0.$$

Full text

spiral wave chaos, alternans, conduction block, cardiac arrhythmia, fibrillation

Atrial fibrillation is a type of cardiac arrhythmia featuring multiple wavelets that continually interact with each other, appear, and disappear. The genesis of this spatiotemporally chaotic state has been linked to the alternans instability that leads to conduction block and wave breakup, generating an increasing number of wavelets. Less clear are the dynamical mechanisms that sustain this state and, in particular, maintain the balance between the creation and destruction of spiral wavelets. Even the relation between wave breakup and conduction block, which is well-understood qualitatively, at present lacks proper quantitative description. This paper introduces a topological description of spiral wave chaos in terms of the dynamics of wavefronts, wavebacks, and point defects – phase singularities – that anchor the wavelets. This description both allows a dramatic simplification of the spatiotemporally chaotic dynamics and enables quantitative prediction of the key properties of excitation patterns.

I Introduction

Atrial fibrillation (AF) is the most common sustained cardiac arrhythmia Waktare and Camm (1998). While not itself lethal, it has a number of serious side effects, such as increased risk of stroke and systemic thromboembolism Noel et al. (1991). The origin of AF has been debated through much of the previous century Nattel, Li, and Yue (2000). In 1913, Mines proposed that fibrillation is caused by a reentrant process Mines (1913), which leads to a high-frequency wave propagating away from the reentry site and breaking up into smaller fragments. This mechanism is presently referred to as anatomical reentry and requires a structural heterogeneity of the cardiac tissue, such as a blood vessel (e.g., vena cava). Reentry could also be functional Allessie, Bonke, and Schopman (1977), where the heterogeneity (i.e., tissue refractoriness) is dynamical in nature. Neither picture, however, explains the complexity and irregularity of the resulting dynamics.

The first qualitative explanation of this complexity came in the form of the multiple wavelet hypothesis proposed by Moe Moe (1962). In this hypothesis multiple independent wavelets circulate around functionally refractory tissue, with some wavelets running into regions of reduced excitability and disappearing and others breaking up into several daughter wavelets, leading to a dynamical equilibrium. This picture was subsequently evaluated and refined based on numerical simulations Moe, Rheinboldt, and Abildskov (1964) and experiments Allessie (1985).

Krinsky Krinskii (1966) and then Winfree Winfree (1987) suggested that the dynamical mechanism of fibrillation relies on the formation and interaction of spiral waves. Spiral waves rotate around phase singularities that may or may not move, producing reentry that requires neither structural nor dynamical heterogeneity in refractoriness. The presence and crucial role of spiral waves in AF was later confirmed in optical phase mapping experiments Davidenko et al. (1992); Pertsov et al. (1993); Chen et al. (2000). Experimental evidence shows that spiral waves tend to be very unstable: only a small fraction of these complete a full rotation Chen et al. (1998); Lee et al. (2001); Liu et al. (2003) with the dynamics dominated by what appears to be wavebreaks or wave breakups (WBs).

Although there is plentiful experimental and computational Fenton et al. (2002) evidence that WBs play a crucial role in the transition to fibrillation, it is far less obvious that this mechanism is essential for maintaining AF. As Liu et al. Liu et al. (2003) write, “for a wave to break, its wavelength must become zero at a discrete point somewhere along the wave. This can happen if the wave encounters refractoriness that creates local block (wavelength = 0), while propagating elsewhere. Therefore, WBs can be detected at locations where activating wavefronts meet the repolarization wavebacks.” The data produced by experimental studies is highly unreliable in this regard, since detecting the position of wavefronts and wavebacks based on optical recordings is far from straightforward. Numerical simulations, on the other hand, have focused mostly on the transition, rather than sustained AF. Theoretical studies of model systems such as the complex Ginzburg-Landau equation Aranson and Kramer (2002) and FitzHugh-Nagumo equation Panfilov and Hogeweg (1993) lack dynamical features, such as the alternans instability Nolasco and Dahlen (1968), that are believed to play an essential role in conduction block that leads to WBs and AF Narayan et al. (2011).

Even if WBs do play a role in maintaining AF, they tell only a part of the story. Indeed, in sustained AF, despite some variation, the quantitative metrics such as the number of wavelets or phase singularities, have to remain in dynamical equilibrium. While WBs may explain the increase in the number of wavelets and phase singularities, it cannot explain how these numbers might ever decrease. The multiple wavelet hypothesis Moe (1962) comes the closest to providing all the necessary ingredients for such a dynamical equilibrium, but it lacks sufficient detail to be either validated or refuted.

The main objective of this paper is to construct a mathematically rigorous topological description of the dynamics and to use this description to characterize and classify different dynamical events that change the topological structure of the pattern of excitation waves in a state of sustained atrial fibrillation. We will focus on the smoothed version Byrne, Marcotte, and Grigoriev (2015) of the Karma model Karma (1993, 1994),

[TABLE]

where $\mathbf{u}(t,\mathbf{x})=[u_{1},u_{2}](t,\mathbf{x})$,

[TABLE]

where $\Theta_{s}(u)=[1+\tanh(su)]/2$, $u_{1}$ is the (fast) voltage variable, and $u_{2}$ is the (slow) gating variable. The parameter $\epsilon$ describes the ratio of the corresponding time scales, $s$ is the smoothing parameter, and the diagonal matrix $D$ of diffusion coefficients describes spatial coupling between neighboring cardiac cells (cardiomyocytes). The parameters of the model Marcotte and Grigoriev (2016a) are $M=4$, $\epsilon=0.01$, $s=1.2571$, $\beta=1.389$, $u_{1}^{*}=1.5415$, $D_{11}=4.0062$, and $D_{22}=0.20031$. This is the simplest model of cardiac tissue that develops sustained spiral wave chaos from an isolated spiral wave through the amplification of the alternans instability resulting in conduction block and wave breaks, mirroring the transition from tachycardia to fibrillation.

The outline of the paper is as follows: Section II introduces the topological description of the complicated multi-spiral states. Section III discusses the relationship between tissue refractoriness and conduction block. The dynamical mechanisms that maintain spiral wave chaos are presented in Section IV. Section V discusses the statistical measures quantifying sustained dynamics, and Section VI contains the discussion of our results and conclusions.

II Wave anatomy

To quantitatively describe the topological changes such as wave breakup and creation/destruction of spiral cores we must first discuss the anatomy of excitation waves and define the appropriate terminology.

II.1 Wavefront and waveback

The region of excitation can be thought of as being bounded by the wavefront, which describes fast depolarization of cardiac cells, and waveback, which describes a typically slower repolarization. The most conventional definition of the wavefront and waveback that has been used in both experiment Lee et al. (2001); Roberts, Hersh, and Scher (1979) and numerical simulations Fenton et al. (2002); Murthy et al. (2013) is based on a level-set of the voltage variable, $u_{1}(t,\mathbf{x})=\bar{u}_{1}$. If we define the region of repolarization

[TABLE]

then the wavefront/back is the part of the level set outside/inside $R$.

The choice of the voltage threshold $\bar{u}_{1}$ is arbitrary and is typically taken as a percentage of the difference between the voltage maximum and its value in the rest state Fenton et al. (2002). This very simple definition allows easy identification of the action potential duration (APD) and diastolic interval ($\mathrm{DI}$), where the percentage is often used as a subscript which refers to the choice of the threshold (e.g., $\mathrm{APD}_{80}$ corresponds to 80% of the difference).

Analytical studies tend to use a different definition Tyson and Fife (1980); Tyson and Keener (1988); Cross and Hohenberg (1993) for the wavefront and waveback which is based on scale separation between the dynamics of fast variables, such as voltage, and slow variables, such as potassium concentration (the voltage $u_{1}$ and the gating variable $u_{2}$, respectively, in the Karma model). For the simplest two-variable models (Karma, Barkley Barkley (1991), FitzHugh-Nagumo FitzHugh (466), Rinzel-Keller Rinzel and Keller (1973), etc.), in the limit $\epsilon\to 0$, the excitation waves are related to a limit cycle oscillation in the $(u_{1},u_{2})$ plane governed by the system of coupled ordinary differential equations Cross and Hohenberg (1993)

[TABLE]

The wavefront and waveback correspond to the segments of the limit cycle solution of (II.1) connecting the two stable branches of the $u_{1}$-nullcline $f_{1}({\bf u})=0$ for which $\mathrm{d}u_{2}/\mathrm{d}u_{1}<0$ (they are denoted with superscripts $-$ and $+$, as illustrated in Fig. 1(a)).

In the limit $\epsilon\to 0$ these segments become horizontal lines and describe very fast (in time) variation of the voltage variable. In space, both the wavefront and the waveback have widths that scale as $\sqrt{\epsilon}$. In the limit $\epsilon\to 0$ they become very sharp (green curves in Fig. 1(b)) and can be thought of as the boundaries of the region of excited tissue

[TABLE]

shown as the green shaded area in Fig. 1(b). These boundaries can be defined with equal precision using any curve in the $(u_{1},u_{2})$ plane that bisects both the wavefront and the waveback segments of the limit cycle. If we define this curve as the zero level set of an indicator function

[TABLE]

the wavefront and the waveback in the physical space at a particular time $t$ are given by

[TABLE]

In particular, a level set of the voltage variable discussed previously corresponds to a vertical line in Fig. 1(a), $g(\mathbf{u})=u_{1}-\bar{u}_{1}$.

A less arbitrary and dynamically better justified choice $g(\mathbf{u})=f_{1}^{0}({\bf u})$ corresponds to the unstable branch of the $u_{1}$-nullcline (for which $\mathrm{d}u_{2}/\mathrm{d}u_{1}>0$). In order to generalize this choice to finite values of $\epsilon$, the unstable branch has to be extended beyond its end points ${\bf u}^{\lessgtr}$ where $\mathrm{d}u_{2}/\mathrm{d}u_{1}=0$, e.g.,

[TABLE]

The excited region where $g({\bf u})>0$ according to (8) is shaded green in Fig. 1(a). Yet another alternative suggested by Fig. 1(a) is to use the $u_{2}$-nullcline $g(\mathbf{u})=f_{2}({\bf u})$ that also bisects both the wavefront and the waveback.

As time evolves, the level set (7) moves with normal velocity

[TABLE]

where $\mathbf{b}=\partial g/\partial\mathbf{u}$, $\partial_{n}=\mathbf{n}\cdot\nabla$ is the directional derivative, and $\mathbf{n}$ is the outside normal to $\partial E$. This normal velocity is taken to be positive for the wavefront and negative for the waveback. In particular, for $g(\mathbf{u})=u_{1}-\bar{u}_{1}$, (9) simplifies, yielding $c=-\partial_{t}u_{1}/\partial_{n}u_{1}$, so the sign of $c$ corresponds to the sign of $\partial_{t}u_{1}$ for proper choices of $\bar{u}_{1}$ such that $\partial_{n}u_{1}<0$ over the entire level set.

Whatever the choice of the bisecting curve $g({\bf u})=0$ in the $(u_{1},u_{2})$ plane is, it may not define a continuous curve in the physical space at all times. Indeed, the PDE model (1) does not take the cellular structure of the tissue into account. Instead, a spatial discretization of (1) should be used, so that the field $\mathbf{u}$ becomes a discontinuous function of space. In this case the wavefront and waveback should instead be defined as the boundary $\partial E$ of the region

[TABLE]

rather than the level set (7).

There are several serious problems with the “local” definitions discussed above, which are based on the kinetics of isolated cardiac cells. For one, they ignore the coupling between neighboring cells in tissue (electrotonic effects) and hence cannot correctly describe the essential properties of excitability and refractoriness, making quantitative description of conduction block impossible. Furthermore, since the value of $\epsilon$ is not vanishingly small for physiologically relevant models, the widths of both the wavefront and the waveback become finite as well, so different choices of $g(\mathbf{u})$ can produce rather distinct results in the physical space. Additional complications will be discussed below.

The definitions of the wavefront and waveback can be generalized for a tissue model by noticing that the level sets such as $f_{1}^{0}(\mathbf{u})=0$ or $f_{2}(\mathbf{u})=0$, in the limit $\epsilon\to 0$, coincide with the level sets $\partial_{t}u_{1}=0$ and $\partial_{t}u_{2}=0$. These are special cases of a more general relation

[TABLE]

where $\mathbf{a}=(\cos\alpha,\sin\alpha)$ and $0<\alpha<\pi$ is a parameter that can be chosen to properly describe the refractoriness and excitability of the model for finite values of $\epsilon$. The wavefront and waveback can again be distinguished as the parts of the level set that lie outside or inside $R$, respectively. We will set $\alpha=\pi/2$ below, which yields the following definition

[TABLE]

The level set (12) is shown as the dashed blue line in Fig. 1(b). As Fig. 2(a-b) illustrates, for the Karma model, it gives an extremely good agreement with the more conventional definition based on (8) for both the wavefront and the waveback. Variation in $\alpha$ by $O(\epsilon)$ has a very weak effect on the position of the wavefront (which is very sharp), but has a more pronounced effect on the position of the waveback (which is much broader).

II.2 Phase singularities

Typically (although certainly not always Biktashev (2005)), the temporal frequency and wavelength of spiral waves are controlled by their central region, usually referred to as a spiral core or rotor Winfree (1978). This region is spatially extended and its size can be characterized using the adjoint eigenfunctions of the linearization Biktasheva, Elkin, and Biktashev (1998); Marcotte and Grigoriev (2016b). In practice it is more convenient to deal with a single point that describes the location of the core region. In particular, the center of this region is associated with a phase singularity, where the amplitude of oscillation vanishes. The location of the phase singularity depends on the definition of the phase, however, and the proper definition is far from obvious for strongly nonlinear oscillations characteristic of excitable systems. The methods based on phase Gray, Pertsov, and Jalife (1998); Iyer and Gray (2001) or amplitude Byrne, Marcotte, and Grigoriev (2015) reconstruction rely on the dynamics being nearly recurrent and break down for spatiotemporally chaotic states featuring frequent topologically nontrivial events such as the creation or annihilation of spiral cores.

A more conventional (and convenient) approach is to use instead the spiral tip, which is a point on the boundary $\partial E$ of the excited region that separates the depolarization wavefront from the repolarization waveback (green circle in Fig. 1(b)). A number of different definitions of the spiral tip have been introduced in the literature. The most popular are the ones based on the level-set intersection (LSI) Barkley, Kness, and Tuckerman (1990)

[TABLE]

or zero normal velocity (ZNV) Fenton and Karma (1998)

[TABLE]

or the curvature $\kappa$ of the level set $u_{1}(t,\mathbf{x})=\bar{u}_{1}$ Beaumont et al. (1998). In particular, ZNV and LSI define the spiral tip(s) $\mathbf{x}_{p}(t)$ as intersection of two level sets which are much easier to compute than the curvature. LSI can be thought of as the limiting case of ZNV where $D_{11}\to 0$, so the difference between the positions of spiral tips defined using these two methods provides a measure of the importance of electrotonic effects.

Albeit they could be simple to define, spiral tips typically exhibit spurious dynamical effects. For instance, they move (in circular trajectories) for spiral wave solutions of (1) rigidly rotating around an origin $\mathbf{x}^{\prime}$ (née relative equilibria) which satisfy

[TABLE]

where $\partial_{\theta}=\hat{\mathbf{z}}\cdot(\mathbf{x}-\mathbf{x}^{\prime})\times\nabla$ and $\omega=2\pi/T$ is the angular frequency. Indeed, even if the normal velocity (9) of the spiral tip vanishes, its tangential velocity will not vanish, unless $\bar{u}_{1}=u_{1}(\mathbf{x}^{\prime})$. Hence, spiral tips are not ideally suited to be used as indicators of the wave dynamics.

The phase singularity, unlike the spiral wave tip, should remain stationary for a rigidly rotating spiral wave. This requires that the location $\mathbf{x}_{o}(t)$ of every phase singularity satisfies

[TABLE]

Equivalently, $\mathbf{x}_{o}(t)$ correspond to the intersections of the level sets

[TABLE]

and $\partial E$ defined according to (12), i.e., phase singularities are points on the boundary of the excited region that separate the refractory region from the excitable region. Note that, for $\partial E$ defined by (11), its intersections with $\partial R$ are independent of $\alpha$, and so is the definition of the phase singularities. This is explicit in the definition (16).

It is easy to see that the boundaries $\partial E$ and $\partial R$ merely correspond to different level sets $\varphi=\alpha\pm\pi/2$ and $\varphi=\pm\pi/2$ of the phase field

[TABLE]

so $\mathbf{x}_{o}(t)$ indeed corresponds to a phase singularity of the phase field (18). Since they are defined locally (just like the spiral wave tips defined via LSI and ZNV), the phase singularities can be easily determined for arbitrarily complicated solutions. More generally, $\mathbf{x}_{o}(t)$ can be interpreted as the instantaneous center of rotation for slowly drifting spiral waves, i.e., for which the rotation-averaged translation of the spiral wave core is much smaller than the typical propagation velocity $c$ of excitation waves.

The positions of spiral wave tips and phase singularities are compared in Fig. 2(c). For the LSI and ZNV definitions, the positions of the spiral tips are shown for $1.68\leq\bar{u}_{1}\leq 2.11$, corresponding to the voltage threshold between $\mathrm{APD}_{70}$ and $\mathrm{APD}_{90}$, respectively Fenton et al. (2002). Clearly, electrotonic effects are non-negligible for the present model, as the tip positions predicted by LSI are dramatically different from those defined by ZNV over a range of choices of $\bar{u}_{1}$. The tip positions defined by ZNV much more closely match the phase singularities defined by (16), and mostly differ in position in the direction normal to the wavefront (along the local gradient of $u_{1}$).

In conclusion of this section, we should mention that, for a typical multi-spiral solution, $u_{1}$ varies rather significantly across the spatial domain, while $u_{2}$ is restricted to a narrow range of $O(\epsilon)$ width around the value $\bar{u}_{2}$ that corresponds to the “stall solution” for a planar front connecting the two stable branches of the $u_{1}$-nullcline

[TABLE]

where $f_{1}^{\pm}(u_{1}^{\pm},\bar{u}_{2})=0$ (cf. Fig. 1(a)). This is a characteristic value that corresponds to the phase singularities, as shown by Fife Fife (1984). For the parameters used here Marcotte and Grigoriev (2016b) $\bar{u}_{2}=0.9724$.

II.3 Topological description

We can associate chirality (topological charge) $q_{j}=\pm 1$ with each of the phase singularities (enumerated by $j=1,2,\dots$), which determines whether the spiral wave rotation is counter- or clockwise. Chirality can be defined locally Li, Deng, and Zhang (2013) as

[TABLE]

For spatially discretized models it is more reliable to use a nonlocal definition of chirality. Let us define a neighborhood of each phase singularity $\mathbf{x}_{o,j}$ using the window function

[TABLE]

where $r_{j}=|\mathbf{x}-\mathbf{x}_{o,j}|$ and $d_{j}=\min_{k\neq j}|\mathbf{x}_{o,j}-\mathbf{x}_{o,k}|$ is the distance to the nearest distinct phase singularity. Further, let us define the pseudo-chirality $\tilde{q}_{j}$ of each spiral wave as the value for which the function

[TABLE]

is minimized. The functional $J(\tilde{q}_{j})$ defines a local reference frame rotating with angular velocity $\tilde{q}_{j}\omega$ around the phase singularity $\mathbf{x}_{o,j}$; this functional is minimized for spiral waves that are stationary in that reference frame. For a single rigidly rotating spiral wave, chirality is precisely $\pm 1$. Minimization of (22) for complex multi-spiral states produces pseudo chirality values equal to $\pm 1$ within a few percent, such that we can safely define $q_{j}={\rm sign}(\tilde{q}_{j})$. In practice, this definition proves very robust when spiral cores are sufficiently well separated, i.e., when $d_{j}$ exceeds the width of an isolated spiral core Marcotte and Grigoriev (2016a).

Since phase singularities by definition (16) lie on the level set $\partial E$, for periodic boundary conditions, wavefronts and wavebacks can only terminate at a spiral core (or, more precisely, phase singularity). Conversely, in multi-spiral states, each wavefront and waveback is bounded by a pair of spiral cores of opposite chirality. In modern electrophysiology literature wavelets are identified with wavefronts Chen et al. (2000). Consequently, the events when a wavelet is created (destroyed) are associated with an increase (decrease) in the number of spiral cores by two. Although the total number of spiral cores is not conserved, the total topological charge

[TABLE]

is conserved Davidsen, Glass, and Kapral (2004); Zhang et al. (2007).

A number of topologically distinct processes which respect (23) are possible. Although some of these correspond to the time-reversed version of the others, the dynamics of the dissipative systems are not time-reversible and do not have to respect the symmetry between these processes. In fact, as we will see below, in excitable systems such as the Karma model the dominant topological processes increasing/decreasing the number of spiral cores are not related by time-reversal symmetry.

III Conduction block and tissue refractoriness

As we mentioned previously, conduction block plays a major role in wave breakup, which is essential for transition to fibrillation and spiral wave chaos in general. The origin of conduction block can be structural, i.e., related to tissue heterogeneity Courtmanche and Winfree (1991), but it can also be dynamical, i.e., occur in homogeneous tissue as a result of an instability. For instance, conduction block can occur when a receding waveback is moving slower than the subsequent advancing wavefront Courtemanche (1996). In fact, there is a variety of other dynamical mechanisms leading to conduction block Fenton et al. (2002).

Conduction block refers to the failure of an excitation front to propagate because the tissue ahead of it is refractory and cannot be excited. Refractoriness is traditionally Merideth et al. (1968) defined on the level of individual cardiac cells by quantifying whether a voltage perturbation applied to the quiescent state of the cell will trigger an action potential. These definitions are not particularly useful for understanding conduction block in tissue for two reasons Starmer et al. (1993); Starobin, Zilberter, and Starmer (1994); Neu, Preissig, and Krassowska (1997): First of all, in tissue the excitation wave is triggered by coupling between neighboring cells, rather than a voltage perturbation to an isolated cell. Second, in tissue, especially during AF or tachycardia, which are both characterized by very short $\mathrm{DI}$s, the cells never have sufficient time to return to the rest state.

III.1 Low-curvature wavefront

In the Karma model, conduction block can arise as a result of discordant alternans instability Echebarria and Karma (2002, 2007) which leads to variation in the width and duration of action potentials. For excitation waves with low curvature, we can determine the boundary of the refractory region by considering a one-dimensional periodic pulse train. In the reference frame moving with velocity $c$ of the wavefront, the voltage variable satisfies the evolution equation

[TABLE]

provided the pulse train does not change shape, where $u^{\prime}_{1}=\partial_{\xi}u_{1}(\xi)$, $u^{\prime\prime}_{1}=\partial_{\xi}^{2}u_{1}(\xi)$, and $\xi=x-ct$. For sufficiently small $\mathrm{DI}$, the conduction velocity $c$ decreases monotonically with $\mathrm{DI}$ and vanishes identically at finite $\mathrm{DI}$ Karma (1994). This means that there are no propagating solutions below this value of $\mathrm{DI}$. At the critical value of $\mathrm{DI}$, we have $c=0$, so the wavefront fails to propagate when $D_{11}u^{\prime\prime}_{1}+f_{1}(\mathbf{u})=0$. For plane waves in two dimensions, $u^{\prime\prime}_{1}=\nabla^{2}u_{1}$, so combining this with the evolution equation (1) we find that the boundary of the refractory region is given by

[TABLE]

and coincides with the boundary (17) of the repolarization region. Similarly, the refractory region can be identified with the region of repolarization (3). This makes intuitive sense: whatever the conditions are, the voltage increases outside the refractory region.

Although derived for a very special case of one-dimensional periodic pulse trains, this definition of the refractory region works well even for states that are not time-periodic and feature excitation waves with significant curvature. This is illustrated in Fig. 3 which shows the time trace of of the variables $u_{1}(t,\mathbf{x}_{0})$ and $u_{2}(t,\mathbf{x}_{0})$ for a spatiotemporally chaotic solution similar to that shown in Fig. 2. The point $\mathbf{x}_{0}$ was chosen near the spatial location where conduction block occurs (such as the center of the marked region in Fig. 7 below). The excited and refractory intervals (temporal analogues of the excitable and refractory regions) are shown as red- and blue-shaded rectangles in Fig. 3; they are bounded by the level sets $\partial R$ and $\partial E$ in Fig. 7.

As expected, we find that, when the wavefronts are well separated from the trailing edges of the refractory intervals (e.g., at $t\approx 140$), long, large-amplitude action potentials are found. This is in sharp contrast with the short and low-amplitude action potential that is initiated at $t\approx 29.5$, soon after the previous refractory interval ends at $t\approx 26.5$. As the location of conduction block is approached (not shown), the separation between the trailing edge of the refractory interval and the subsequent excitation wavefront vanishes along with the action potential itself. This suggests that conduction block occurs when and where the level sets $\partial R$ and $\partial E$ first touch, in agreement with Winfree’s critical point hypothesis Winfree (1989).

III.2 High-curvature wavefront

There are, however, other dynamical mechanisms that can lead to conduction block. Consider, for instance, the opposite situation when the curvature of the wavefront is high. For curved wavefronts the propagation speed $c$ decreases as the curvature $\kappa$ of the wavefront increases, so there is a critical value of the curvature at which the wave fails to propagate Fast and Kléber (1997). It can be estimated in the limit $D_{22}\rightarrow 0$ using the eikonal approximation Keener (1986) which gives

[TABLE]

where $c_{0}$ is the velocity of a planar wavefront. Using the value of $c_{0}=\lambda/T\approx 1.44$ which corresponds to a large rigidly rotating spiral wave Marcotte and Grigoriev (2015), we find $\kappa^{-1}=r_{c}\approx 2.8$, where $r_{c}$ is the critical radius of curvature of the wavefront.

A more accurate estimate for $r_{c}$ can be obtained using the condition (25) for conduction block

[TABLE]

and the definition of the wavefront (12) rewritten via (1)

[TABLE]

in polar coordinates $(r,\theta)$. For small $\epsilon$, $u_{2}$ varies slowly in both time and space and can be considered a constant, $u_{2}=\bar{u}_{2}$ given by (19). The first three terms of (28) can therefore be ignored, so (28) reduces to $f_{2}(\mathbf{u})=0$. The third term $D_{11}r^{-2}\partial_{\theta}^{2}u_{1}$ in (27) can also be neglected, since $u_{1}$ varies much faster in the direction normal to the wavefront $r=r_{c}$ than in the tangential direction. Since $u_{2}$ is constant, subject to boundary conditions $\partial_{r}u_{1}=0$ at $r=0$ and $r=\infty$, (27) has rotationally symmetric solutions $u_{1}=u_{1}(r)$ with a stationary wavefront at the critical radius $r_{c}$ given by

[TABLE]

The stationary solution of (27) and (29) is shown in Fig. 4. It corresponds to the critical radius of curvature $r_{c}\approx 6$, which is a factor of two larger than the value obtained using the eikonal approximation (26). This solution is a two-dimensional analogue of the one-dimensional “critical nucleus” for excitation Neu, Preissig, and Krassowska (1997). Wavefronts with radius of curvature larger than $r_{c}$ propagate forward, while the wavefronts with radius of curvature smaller than $r_{c}$ retract (i.e., become wavebacks).

Before we discuss the numerical results, let us emphasize that, with the proper choice of variables, the definitions of the wavefronts and wavebacks (12), leading and trailing edges of the refractory region (17), and phase singularities (16) are model-independent and can be used to analyze both numerical and experimental data, provided that measurements of two independent variables (e.g., voltage and calcium) are available. Generalization of the topological description presented above to higher-dimensional models is discussed in the Appendix.

IV Numerical results

As mentioned in the introduction, during AF most spiral waves do not complete a full rotation. Spiral wave chaos in the Karma model produces qualitatively similar dynamics: topological changes involving a change in the number of spiral wave cores occur on the same time scale as the rotation. (For parameters considered in this study rotation period is $T\approx 51$, which corresponds to $127$ ms in dimensional units Marcotte and Grigoriev (2015).) The larger the spatial domain, the more frequent are the topological changes in the structure of the solution. However, as we discussed in Sect. II.3, each topological event is essentially local and involves either birth or annihilation of a pair of spiral cores of opposite chirality.

Of the different types of topological events, spiral wave breakup – associated with a creation of a new pair of spiral cores – received a lion’s share of the attention due to its role in the initiation of fibrillation. However, the number of spiral cores cannot increase forever; eventually a dynamic equilibrium is reached when the number of cores fluctuates about some average, with core creation balanced by core annihilation. To the best of our knowledge, the process(es) responsible for core annihilation, however, have never been studied systematically. To investigate which of the topological events dominate and what the dynamical mechanisms underlying these events are, we performed a numerical study of the Karma model (1)-(2) on a square domain of side-length $L=192$ ($5.03$ cm), which is close to the minimal size required to support spiral wave chaos. Spatial derivatives were evaluated using a second-order finite-difference stencil and a fourth-order Runge-Kutta method was used for time integration Marcotte and Grigoriev (2015). To avoid spurious topological transitions involving a boundary, periodic boundary conditions were used, unless noted otherwise.

Before identifying topological transitions in the numerical simulations, it is worth enumerating the topologically distinct local configurations. In the following it will be convenient to use the following shorthand notations: $\partial E^{+}$ (wavefront, $\partial_{t}u_{1}>0$), $\partial E^{-}$ (waveback, $\partial_{t}u_{1}<0$), $\partial R^{+}$ (leading edge of the refractory region, $\partial_{t}u_{2}>0$), and $\partial R^{-}$ (trailing edge of the refractory region, $\partial_{t}u_{2}<0$). For a wave train in the region with no spiral cores, the boundaries of the refractory and excitable regions will follow the periodic sequence $(\dots$, $\partial R^{-}$, $\partial E^{-}$, $\partial R^{+}$, $\partial E^{+}$, $\dots)$ in the co-moving frame (cf. Fig. 5(e)).

IV.1 Virtual pairs

Deformation of (nearly planar) waves due to instability or heterogeneous refractoriness can lead to intersection of any pair of adjacent level sets (e.g., $\partial E^{+}$ and $\partial R^{-}$) and, correspondingly, creation of a new pair of spiral cores. From the topological perspective, there are four distinct possibilities shown in Figs. 5(b), 5(d), 5(f), or 5(h). The corresponding configurations are all transient and only persist for a fraction of the revolution time $T$ of a typical spiral wave during topological transitions. Each of these transient configurations can undergo a total of five distinct topological transitions, including a reverse transition back to the initial configuration with nonintersecting level sets, associated with the destruction of the two new phase singularities. The other four persistent possibilities will be considered in subsequent sections.

The transitions that correspond to crossing of two level sets followed by the reverse transition (indicated by horizontal or vertical double gray arrows in Fig. 5) produce a “virtual” core pair that appears and quickly disappears, restoring the original topological structure. The number of spiral cores and wavelets before and after these transitions remains exactly the same, so while such events do occur rather frequently they do not play a dynamically important role and can be safely ignored. As discussed below, the topological transitions identified with the arrows in Fig. 5 are all very fast; they occur on a time scale much shorter than the typical rotation period of a spiral. The much slower transitions between panels 5(a) $\to$ 5(b) $\to$ 5(c) $\to$ 5(f) $\to$ 5(i) $\to$ 5(h) $\to$ 5(g) $\to$ 5(d) $\to$ $\cdots$ associated with figure-8 re-entry for two well-separated counter-rotating spirals are not shown for clarity. While the topology of some level sets (either $\partial E$ of $\partial R$) changes during these transitions, neither the topological charge nor the number of phase singularities does, so these are not proper topological transitions, as defined in Section II.3. In contrast, the topology of $\partial E$ and $\partial R$ may not change during the proper topological transitions.

The trajectories of the two phase singularities and the distance between them in a representative example of a virtual pair are shown in Fig. 6. The phase singularities do not move far from their initial positions and remain close at all times. In fact, the distance between them never exceeds a fraction of the typical separation $d_{0}$ between persistent spiral cores (to be be discussed in more detail in Sect. V). Since the level sets are smooth curves, their intersections that define the positions of the cores move with infinite velocities at the time instants when the cores are created and destroyed. As a result, their motion near those times can not be resolved using time stepping, explaining the gaps at the beginning and the end of the trajectories in Fig. 6(a).

The transient configuration shown in Fig. 5(f) deserves a special mention. It corresponds to the phenomenon of “back-ignition” observed in some reaction-diffusion models whereby the waveback can become a source of a new backward propagating wave under appropriate conditions. While topologically permissible, this configuration is only observed as a very short transient in the present model, reflecting the asymmetry of initial conditions imposed by the dynamics. The relative likelihood of the four transient configurations shown in Figs. 5(b), 5(d), 5(f), and 5(h) can in principle be computed using stability analysis of a planar wave train solution for any tissue model, but this is outside the scope of the present study.

IV.2 Wavelet/pair creation

Next consider the transitions from the four intermediate configurations shown in Figs. 5(b), 5(d), 5(f), and 5(h) to configurations other than the initial one shown in Fig. 5(e). For each of these four transient configurations there are four distinct possibilities. Two possibilities are shown in Fig. 5 (the other two will be discussed in the next Section): either of the two level set fragments connecting the cores can reconnect with the neighboring level set of the same type. For instance, the configuration shown in Fig. 5(h) can transform to the configuration shown in Figs. 5(g) or 5(i). If the crossing and reconnection occur simultaneously, the transition occurs directly between the persistent configuration shown in Fig. 5(e) and one of the persistent configurations shown in Figs. 5(a), 5(c), 5(g) or 5(i) without passing through any of the intermediate transient configurations. The dynamically allowed direct transitions, as determined based on the results of numerical simulations, are shown as diagonal gray arrows.

Note that the transition between the configurations shown in Fig. 5(e) and 5(a) corresponds to wave breakup. It occurs when and where the wavefront reconnects with the waveback of the same excitation wave Fife (1985); Winfree (1989) as a result of conduction block. Since this topological process increases the number of disconnected excited regions, it is quite natural to find that it plays an important role in the transition from, say, normal rhythm or tachycardia (featuring a single excitation wave) to AF (featuring many separate wavelets). While wave breakup may be prevalent during the initial stage when AF is being established, we have not found a single instance of this topological event in our numerical simulations of sustained spiral wave chaos, casting serious doubt on the premise that wave breakup plays a dynamically important role in maintaining AF in tissue or in other models.

Our numerical simulations reveal only one topological process that leads to a lasting increase in the complexity of the pattern. This process which we call “wave coalescence” corresponds to the transition from the initial configuration shown in Fig. 5(e) to the final configuration shown in Fig. 5(i), either directly or through the intermediate transient configurations shown in Figs. 5(f) and 5(h). A representative example from the simulations is shown in Fig. 7, where a purple rectangle marks the region of interest. Outside of this region the wavefronts are well-separated from the refractory tails, but inside the separation is markedly smaller (cf. Fig. 7(a)). The separation quickly decreases (cf. Fig. 7(b)) until the level sets $\partial E^{-}$ and $\partial R^{-}$ cross and two new spiral cores with opposite chirality are created (cf. Fig. 7(c)). Immediately after this the two parts of the level set $\partial E^{+}$ reconnect, bringing the configuration to the topological state shown in Fig. 5(i). The cores separate (cf. Fig.  7(c)) and the excited regions of two subsequent waves coalesce in the gap flanked by these two cores.

Due to the high curvature of $\partial E$ the two new cores are quickly pulled apart, and two new counter-rotating spiral waves emerge, “locking in” the resulting topological configuration. This is illustrated by Fig. 8 which shows the trajectories of the cores and the distance between them. It is worth noting that, before the spiral waves complete even half a revolution, the separation between the cores approaches the typical equilibrium distance Marcotte and Grigoriev (2016b) $d_{0}$.

Our numerical simulations did not produce any examples of topological transitions to the configurations shown in Figs. 5(c) or 5(g). In the horizontal band bounded by the cores (indicated by lighter-shade gray), the corresponding states are characterized by the voltage variable that is changing slowly in space, since the distance $d$ between the minimum of $u_{1}$ (dashed white line $\partial R^{-}$) and the maximum of $u_{1}$ (solid white line $\partial R^{+}$) is extremely large. Hence, the term $D_{11}\nabla^{2}u_{1}\propto d^{-2}$ in (1) is negligible. Since $D_{22}$ is small, we can ignore the diffusive terms $D_{22}\nabla^{2}u_{2}$ as well and consider all cells in this region to be spatially decoupled, such that their dynamics is described well by (II.1) and the phase diagram shown in Fig. 1(a). Consider the part of the band where $u_{1}$ is slowly and monotonically increasing in time (to the left of $\partial R^{+}$ in Fig. 5(c)) or decreasing in time (to the left of $\partial R^{-}$ in Fig. 5(g)). The cells in this region should be in the state that lies close to either one of the stable $u_{1}$-nullclines ($f_{1}^{-}=0$ in the former case and $f_{1}^{+}=0$ in the latter case). According to Fig. 1(a), this entire region should lie either to the left or to the right of the $u_{2}$-nullcline, so $\partial_{t}u_{2}$ should be sign-definite, while in both Fig. 5(c) and 5(g) the sign of $\partial_{t}u_{2}$ changes (when the level set $\partial E$ connecting the two phase singularities is crossed). Hence, while these configurations are not forbidden on topological grounds, they are forbidden dynamically.

Furthermore, we have not observed transitions from the persistent configurations shown in Figs. 5(a), 5(c), 5(g), and 5(i) to either the transient configurations shown in Figs. 5(b), 5(d), 5(f), and 5(h) or the persistent configuration shown in Fig. 5(e). While these transitions are allowed topologically, they appear to be forbidden dynamically. The only dynamically allowed (direct or indirect) transition irreversibly transforms the configuration with no spiral cores (Fig. 5(e)) to the configuration (Fig. 5(i)) with two spiral cores, increasing the total number of cores by two and the number of wavelets by one.

IV.3 Wavelet/pair destruction

Finally, let us consider the topological transitions from the four intermediate configurations that have not been considered in the previous sections. We have redrawn these four configurations in Fig. 9 in the same locations as in Fig. 5, dropping all non-essential level sets. For each of these intermediate configurations there are two possibilities which involve reconnection between the two extended branches of a level set that terminate at the cores. For instance, the configuration shown in Fig. 9(h) can transform to the configurations shown in Figs. 9(g) or 9(i). None of these transitions (shown as horizontal or vertical gray lines), in either the forward or the reverse direction, have been observed in numerical simulations, however.

As we discussed previously, if the crossing and reconnection of the level sets occur simultaneously, the configuration transitions directly between the persistent configuration with no level set intersections (Fig. 9(e)) and one of the persistent configurations with a pair of intersections shown in Figs. 9(a), 9(c), 9(g), and 9(i) without passing through any of the intermediate configurations. The dynamically allowed direct transitions observed in the simulations are shown as diagonal gray arrows. Note that again there is no time-reversal symmetry: only the transitions that destroy the existing core pairs are dynamically allowed. Therefore the observed direct transitions shown in Fig. 9 reduce the net number of spiral cores and wavelets balancing the increase due to wave coalescence.

The configurations shown in Figs. 9(a), 9(c), 9(g), and 9(i) all describe a pair of counter-rotating spiral waves. In particular, the configurations in Figs. 9(c) and 9(g) correspond to multi-spiral states (wavefronts and/or wavebacks connect spiral cores inside and outside the region shown) and hence are quite typical. On the other hand, the configurations in Figs. 9(a) and 9(i) correspond to configurations with a single pair of spirals (wavefronts and wavebacks connect the phase singularities inside the region shown) and are never observed during sustained spiral wave chaos. Consequently, only transitions from the configurations in Figs. 9(c) and 9(g) are found in the simulations, with the vast majority of transitions involving the former configuration.

To understand why and when this transition happens, consider the interaction between a pair of isolated counter-rotating spiral waves separated by distance $d$. (The interaction is short-range, so the presence of other, remote, spiral wave cores does not change the outcome.) Using the approximate mirror symmetry of the configuration, the dynamics can be understood by considering a single spiral interacting with a planar no-flux boundary at a distance $\zeta=d/2$. As we showed previously Marcotte and Grigoriev (2015, 2016b), at large separations the spiral cores can be considered essentially non-interacting, while at smaller separations the equilibrium distance $d$ becomes quantized, with the smallest stable separation Marcotte and Grigoriev (2016b) equal to $d_{0}=2\zeta_{0}\approx 40$ (10.4 mm in dimensional units) for the values of the parameters considered in this study. For separations below some critical distance $d_{c}<d_{0}$, the cores attract each other, eventually colliding and destroying both spiral waves. As the cores approach each other, the wavefront confined between them collapses, so we will refer to this process as “wavelet collapse” or “wave collapse”.

The details of wave collapse depend on the relation between the initial phase of the spiral waves and separation between their cores. A very typical example of wave collapse is shown in Fig. 10. In this particular example we find that the curvature of the wavefront becomes quite large before collapse. The curvature at which this happens can be related to the mechanism of conduction block discussed in Sect. III.2. Since the cores are moving relatively slowly prior to wave collapse (cf. Fig. 10(a)), as the wavefront propagates its curvature gradually increases (cf. Fig. 10(b)). The largest value of the curvature is related to the distance between the cores, $\kappa^{-1}\approx d/2$. Once the curvature becomes comparable to the inverse of the critical radius $r_{c}\approx 6$, the wave stops propagating, the cores slide towards each other along the wavefront and annihilate (cf. Fig. 10(c)), the wavebacks merge, and the wave starts to retract (cf. Fig. 10(d)).

This picture predicts that the minimal distance at which the spiral cores with opposite chirality can persist without annihilating with each other is given by $d_{c}=2r_{c}=12$ (3.1 mm in dimensional units). This value is in good agreement with the critical isthmus width (2.5 mm) found for conduction block in isolated sheets of ventricular epicardial muscle with an expanding geometry Fast and Kléber (1997). Our numerical simulations show that the minimal distance is $d_{c}=16$ (4.2 mm), also close to the predicted value. The core trajectories and the distance between them in the example from Fig. 10 are shown in Fig. 11. The initial distance in this case was $d=18$ (4.7 mm), illustrating that, under appropriate conditions, wave collapse can also occur for cores separations somewhat larger than $d_{c}$ (but still less than $d_{0}$).

The transition between the configurations shown in Figs. 9(g) and 9(e) corresponds to merger between two wavefronts that were originally separated by a waveback. Hence, we shall refer to this topological transition as a “wave merger” event. Wave mergers, however, are extremely rare, so a reduction in the total number of cores and wavelets is due almost entirely to wave collapse events. This is similar to the dynamical asymmetry between the wave breakup and wave coalescence events. Therefore, dynamical equilibrium in sustained spiral wave chaos can be understood, at least in the Karma model, as a balance between wave coalescence and wave collapse.

V Dynamical Equilibrium

Although the topological description itself is not quantitative, it helps identify the key dynamical mechanisms, such as wave coalescence and wave collapse, responsible for maintaining AF. This should, in turn, enable a quantitative description of the dynamics in general and dynamical equilibrium in particular and give the answers to open questions that have been debated for a long time. For instance, it is presently not well understood either what the minimal size of tissue is that can sustain AF or what the minimal number of wavelets is in sustained AF. The leading-circle concept Allessie, Bonke, and Schopman (1977) suggests that the number of wavelets that the atria can accommodate should be related to the wavelength. Moe’s computer model Moe, Rheinboldt, and Abildskov (1964) predicted that between 23 and 40 wavelets are necessary for the maintenance of AF, while Allessie Allessie (1985) places the minimal number of wavelets between four and six.

These hypotheses can be easily tested in the context of the Karma model. Let us start by determining whether the wavelength ($\lambda=78$ for the values of parameters considered here) is a relevant length scale. The size (diameter) of a reentry circle with the perimeter equal to the wavelength is $d=\lambda/\pi\approx 25$ which is larger that the minimal separation $d_{c}$ between persistent spiral wave cores, but considerably smaller than the minimal stable separation $d_{0}$ between the cores.

To show that $d_{0}$ is the relevant length scale, we computed the probability density function $P(d)$ for core-core separation on a square domain of side $L=192$ (this is the smallest domain with no-flux boundary conditions that supports sustained spiral wave chaos). For each time $t$ and each core $j$ we computed the distance $d_{j}$ to the nearest core (cf. Sect. II.3), then averaged over $j$ and $t$. The resulting distribution, for both no-flux and periodic boundary conditions, is shown in Fig. 12. In both cases we find that the distribution $P(d)$ is rather narrow, with the maximum achieved at $d=d_{0}$.

The effect of the boundary conditions on the shape of the distribution is somewhat subtle: on a bi-periodic domain, the probability of large core separations ($d=O(L)$) is decreased compared with the same size domain with no-flux boundary conditions. Effectively, as there must always be a chirally-matched pair on the periodic domain, the furthest these cores may be is $d_{\rm max}=L/\sqrt{2}$, as opposed to an isolated spiral matched with it’s mirror image across the no-flux boundary, which corresponds to maximal distance $d_{\rm max}=\sqrt{2}L$. Thus, on a periodic domain, the maximal accessible distance is precisely $1/2$ the maximal distance available on a no-flux domain of the same size.

The upper bound for the number of spiral cores can be estimated as the ratio of the total area of the domain (i.e., $L^{2}$) to the area of the smallest tiles Byrne, Marcotte, and Grigoriev (2015) supporting one persistent spiral wave (i.e., $d_{0}^{2}$), that is $\bar{n}_{c}<L^{2}/d_{0}^{2}=23$ (in fact, we should have $\bar{n}_{c}\leq 22$ since the net topological charge is zero). In reality the tiles tend not to be squarish and have a larger area on average, giving a lower average number of spiral cores, $\bar{n}_{c}=10$, as the probability distribution function $P(n_{c})$ illustrates (cf. Fig. 13). The number $n_{w}$ of separate wavelets is exactly half of the number $n_{c}$ of cores (on a domain with periodic boundary conditions), so on average $\bar{n}_{w}=\bar{n}_{c}/2\approx 5$, in perfect agreement with the results of Allessie Allessie (1985). The number of cores exhibits considerable fluctuation (between 4 and 16), correspondingly the number of wavelets varies between 2 and 8. The likelihood of these extreme values is, however, rather small (an order of magnitude smaller than that corresponding to the average value).

The observation that the minimal number of wavelets is just two is a sign that the dynamics are on the border of spontaneous collapse of spiral wave chaos (recall that our domain is just large enough to sustain this regime). We should have $P(0)=0$, because once all the spiral cores disappear, so does the mechanism of reentry (at least in our homogeneous model), resulting in a transition to the rest state or, in the presence of pacing, normal rhythm. The smallest number of spiral cores required for reentry (in a domain with periodic boundary conditions) is two, so one could, in principle, expect $P(2)$ to be nonzero. However, as our results show, the mechanism that sustains spiral wave chaos is wave coalescence, which requires at least two wavelets, and therefore at least four spiral cores, to be present.

VI Discussion

This paper presents a general topological approach for studying spiral wave chaos in two-dimensional excitable media. It is illustrated using the Karma model which, in a certain parameter regime, produces dynamics that are remarkably similar to those observed during atrial fibrillation. Therefore, our results could shed new light on this important and complicated phenomenon.

The confusing and often contradictory results regarding the dynamical origins of AF reported in experimental and numerical studies are to some extent due to the complexity of the patterns of excitation. The descriptive language and intuition developed primarily in the context of simple structures – plane or spiral waves – often fail us when applied to states that are topologically complicated and nonstationary. To give a few examples, the mental picture of a spatially localized excitation wave, or wavelet, that is bounded by a wavefront and a waveback falls apart when applied to complex multi-spiral patterns since the boundary of one excited region is often composed of multiple wavefronts and wavebacks, as Figs. 7 and 10 illustrate. As a result, the number of excited regions almost never corresponds to the number of wavelets. Neither is the notion of a spiral wave immediately useful for describing such complicated patterns, which only resemble spiral waves in small neighborhoods of spiral cores. Similarly, a reduction of complicated field configurations to the number and positions of phase singularities is also problematic, both because they appear, move, and disappear for sustained spiral chaos and because identifying them using existing approaches, such as phase mapping Davidenko et al. (1992); Pertsov et al. (1993), is notoriously unreliable when the data is noisy.

This paper aims to rectify some of these difficulties by introducing a topological description that can rigorously and easily identify the dynamically important elements of the excitation patterns – wavefronts, wavebacks, phase singularities, etc. – without modeling assumptions and in a manner that can be implemented in both simulations and experiments. By defining the phase singularities as intersections of level sets of an appropriately defined phase field, this topological description directly connects the dynamics of excitation waves and phase singularities; it can be used not only to quantify and classify the excitation patterns, but also to identify the dynamical mechanisms that lead to qualitative changes in the pattern. In particular, we show that the qualitative changes can be conveniently described and classified based on the dynamics of spiral cores which are created or destroyed in pairs, leading to an increase or decrease in the number of wavelets, with a one-to-one correspondence between the number of cores and wavelets.

The topological description also allowed us to identify the dominant dynamical mechanisms responsible for maintaining AF in a model of atrial tissue. In particular, it allowed us to make a major discovery with implications that, in all likelihood, go far beyond the simple model considered here. We found that wave breakup due to conduction block that is widely believed to be the key mechanism responsible for maintaining AF plays no role whatsoever in sustaining this regime. While wave breakup does play a key role in the transition to AF, it is a dynamically and topologically distinct event – wave coalescence – that is responsible for maintaining AF. Wave coalescence which leads to the increase in the number of spiral cores and wavelets is balanced by wave collapse which decreases the number of spiral cores and wavelets. It is this delicate balance that is responsible for maintaining the complexity of the pattern and of the dynamics and it is this balance that controls whether AF persists or terminates.

Past studies of the dynamical origins and control of AF tended to focus solely on the mechanism(s) that lead to an increase in the complexity of the pattern. Indeed suppressing the processes that generate new spiral cores and new wavelets is one way to terminate or prevent AF. However, enhancing the processes that destroy the spiral cores and wavelets could be just as effective. Therefore, both wave coalescence and wave collapse are attractive targets for electrical, surgical, and pharmacological approaches to the treatment of AF. While this study has not focused on the interaction of excitation waves with no-flux boundaries, the methods and approaches presented here are applicable to this situation as well. Hence topological analysis could be quite helpful for improving treatment of chronic AF using surgical procedures such as ablation that effectively introduce additional boundaries.

In conclusion, we should point out that our results raise new questions regarding the role of conduction block in maintaining AF. While conduction block undoubtedly plays a crucial role in wave collapse, it is not at all clear that it is relevant in wave coalescence. Therefore, quite paradoxically, we find that conduction block plays a more important role in decreasing the complexity of the excitation pattern than in increasing its complexity. Further studies are needed in order to fully understand the dynamical mechanisms behind wave coalescence, wave collapse, and possibly other topologically allowed events important in maintaining AF using more detailed and physiologically accurate models of atrial tissue.

Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. CMMI-1028133. The Tesla K20 GPUs used for this research were donated by the “NVIDIA Corporation” through the academic hardware donation program. CDM gratefully acknowledges the financial support of the EPSRC via grant EP/N014391/1 (UK).

Appendix A Application to more realistic action potentials

The topological approach introduced here is sufficiently general to be extended to much more complicated electrophysiological models. For example, the definition of the leading and trailing edges of the refractory region (17) relies solely on the voltage variable, while the definition (11) of the wavefronts and wavebacks can be trivially generalized to an $m$-variable reaction-diffusion model by choosing the “weighting” vector $\mathbf{a}=[a_{1},\dots,a_{m}]$ to properly represent the physiological role of different variables in triggering the depolarization front. The phase singularities can then be defined again as the intersection of the boundaries $\partial E$ and $\partial R$, although this does not guarantee that (16) will be satisfied for $m>2$.

For illustration, we used such a generalization to identify the wavefronts, wavebacks, and phase singularities in the four-variable minimal model of Bueno-Orovio et al. Bueno-Orovio, Cherry, and Fenton (2008), considering only the contribution from a single slow variable, $\mathbf{a}=[0,0,1,0]$, as a simple approximation. A virtual pair event shown in Fig. 14 provides an illustration that these definitions are equally successful in describing topological changes in a substantially more complex and detailed model, which is capable of producing quantitatively accurate description of excitation patterns in various types of cardiac tissues, given appropriate choices of parameters (in the present example we used the epicardial parameter set).

The only complication arises when the diffusion coefficients for the slow variable(s) vanish identically, since this can lead to subtle artifacts when discontinuities of the kinetics, e.g., in the switching between on and off states, combine with the high spatial gradients near the wavefront. However, this issue is not a fault of the method, but rather a consequence of the unphysical nature of the simplified ionic kinetics, and can be easily rectified by using an appropriately smoothed version of the model kinetics (as we did in Karma model). Furthermore, even without smoothing, one can simply define the wavefront and waveback as the boundary of the excited region $E$ defined using the indicator function $g({\bf u})=\mathbf{a}\cdot\partial_{t}\mathbf{u}$.