Filament tension and phase-locked drift of meandering scroll waves

Rotating scroll waves are self-organising patterns which are found in many oscillating or excitable systems. Here we show that quasi-periodic (meandering) scroll waves, which include the rotors that organise cardiac arrhythmias, exhibit filament tension when averaged over the meander cycle. With strong filament curvature or medium thickness gradients, however, scroll wave dynamics are governed by phase-locked drift instead of filament tension. Our results are validated in computational models of cycloidal meander and a cardiac tissue model with linear core.

The understanding of the excitation patterns exhibited by circular-core spirals in 2D, and scroll waves in 3D, has much benefited from the analysis of their motion in terms of 'phase singularities', i.e. instantaneous rotation centers for the spirals, Fig. 1, and 'filaments', Fig.  2, for the scrolls [19][20][21]. Much of the theory of meandering spiral waves has been focusing on the origin of the meander bifurcation [10,[22][23][24], which produces epi-or hypocyclodial motion of a spiral tip, as shown in Fig. 1a. However, meandering spirals with 'linear' cores, as in Fig.  1b, may be the building blocks of ventricular fibrillation, which motivated recent work to calculate their leading eigenmodes [25][26][27]. In this Letter, we derive equations of motion for biperiodic meandering 2D spirals and 3D scroll waves, without restriction to a particular shape of meander.
In 3D, it has been shown that the filament of a circularcore scroll wave is characterized by its 'tension' (γ 1 ), which depends on the medium parameters: γ 1 < 0 leads to ever-growing filaments [20,28] if the medium is thick enough [29], resulting in a turbulent, fibrillation-like state, while γ 1 > 0 leads to shrinking of scroll rings, so that only filaments connecting opposite medium boundaries persist. Fig. 2 shows similar behaviour for meandering scroll waves. However, the applicability of the concept of filament tension to meandering scrolls has so far been a conjecture rather than fact. In this Letter, we will show when this is indeed true, and when it is not.
Methods. We investigate spiral-shaped solutions to the reaction-diffusion (RD) system in 2 and 3 spatial dimensions under a small spatiotemporal perturbation h: where u is a column matrix of state variables. Eq.
(1) describes both BZ-like chemical systems and models of cardiac tissue, depending on the choice of the diffusion constants P and reaction kinetics F(u). We consider two different kinetics models: Barkley model Fig. 1a; and the Fenton-Karma (FK) cardiac tissue model with guinea pig (GP) parameters [17], where u = [u, v, w] T , P = diag(0.1, 0, 0) mm 2 /ms, Fig. 1b.
Before presenting our formalism, we recall the classical view on meander in terms of the transition of an unperturbed (h = 0) spiral from rigid to biperiodic rotation via an equivariant Hopf bifurcation [10,[22][23][24]31]. Let ω s be the angular velocity of a rigidly rotating spiral. In a frame of reference rotating with ω s , this spiral is a stationary solution. The corresponding linearised problem always has eigenvalues on the imaginary axis, due to the rotational (λ R = 0) and translational (λ T = ±iω s ) Euclidean symmetry of the plane. If under parameter change another pair of eigenvalues crosses the imaginary axis at λ H = ±iΩ, the solution is becomes time-periodic with period T = 2π/Ω in the rotating frame. The progression of the solution along this cycle can be labeled by the meander phase ψ; by definition, in the absence of perturbation, ∂ t ψ = Ω.
To describe meandering spirals without relying on the proximity of the Hopf bifurcation, we note that in the lab frame, a meandering spiral is a relative periodic orbit, meaning that after the time T the solution returns to the same state up to an orientation-preserving isometry M of the plane [10,23,32,33]. Except for the resonant case, which falls outside our present scope, M is rotation by an angle χ around a point C. After iterating M, the point C emerges as the centre of the meander pattern, see Fig. 1. By construction, χ is defined up to an integer number of full rotations. In terms of the classical approach, we can write χ = ω s T + 2πn, for n ∈ N. For our formalism, the exact choice of χ is not of principal importance. We find it convenient to demand |χ| < π, and define χ as the (smallest) angle between consecutive 'petals' of the tip path, see of reference rotating around C with ω = χ/T , in which the solution is also T -periodic. As before, we define Ω = 2π/T and ∂ t ψ = Ω. Note that this formalism equally holds for both cases shown in Fig. 1.
Let angle φ(t) characterize the orientation of the steadily rotating frame (Fig. 1); by definition ∂ t φ = ω for the unperturbed spiral. The transformation between the lab frame coordinates x a and the rotating frame coordinates ρ A is then where is the rotation matrix over an angle φ. In the rotating frame, Eq. (1) becomes Here τ is time in that frame, ∆ is the Laplacian in the (ρ 1 , ρ 2 ) plane and θ is the polar angle in it, i.e.
The unperturbed meandering spiral wave solution u 0 (ρ 1 , ρ 2 , ψ) to Eq. (2) is 2π-periodic in ψ and satisfies In what follows, we perform a standard perturbation technique used before to derive drift laws for circularcore spiral and scroll waves [19,20,34,35]. This involves linearization of Eq. (3) on u 0 , after which the drift caused by a perturbation h will be given by its projection onto the symmetry eigenmodes.
The linear operatorL associated with Eq. (3) iŝ The operatorL is the same as used for the circular-core case [19,20,35,36]. By differentiating Eq. (3) with respect to ρ 1 , ρ 2 , θ and ψ, we find the four critical eigenmodes: Modes V ± correspond to translations, V φ to rotations and V ψ to shifts in time.
To perform projections onto the symmetry modes, we use the inner products The projectors W M , also called 'response functions (RFs)' in this context [37], are the critical eigenfunctions of the adjoint operatorL † , They are 2π-periodic in ψ, and can be normalized as The biorthogonality property (7) is not practical since it involves averaging over ψ. However, although generally an inner product f | g depends on ψ, for products of eigenfunctions ofL the following Meander Lemma holds (see also [26,33]): If Since this is 2π-periodic, i(λ N − λ M )/Ω needs to be integer. However, in the non-resonant case we have 0 < |χ| < π, whence 0 < |ω/Ω| < 1/2. For the critical modes, (λ N − λ M )/Ω thus cannot be integer, and therefore A M N = 0 and W M | V N = 0. As a corollary, the instant orthogonality (8) also holds for the Cartesian basis of eigenfunctions, i.e. M , N ∈ {ρ 1 , ρ 2 , φ, ψ} and pairs of a critical and non-critical mode.
Results. We are now ready to calculate how a small perturbation h, say of order η, induces spiral wave drift. Still in 2D, we decompose the exact solution to (2) as is made unique at all times by the condition W M |ũ = 0 [35]. Then, we let the frame move with yet unknown perturbed velocities: Finally, projection onto the RFs delivers: The equation of motion (13) describes the spatial drift of the position X 1 , X 2 of the centre of the meander pattern, its orientation φ in the plane and the meander phase ψ of the spiral. It is a fundamental result in this work, as it captures the generic drift response of a meandering spiral wave to small external disturbances h( r, t). Its form was stated before based on symmetry for particular cases of h [24,32,38,39], but without the overlap integral which is necessary to quantitatively predict spiral wave drift.
In this work we choose to further study which has several applications. E.g. in chemical systems, u is a vector of concentrations of reagents, and Eqs. (1), (14) may describe the 'electrophoretic' drift of spiral waves in a constant electrical field E, if Q is the diagonal matrix of electrical mobilities of the reagents. More generically, in any RD system describing 3D scroll waves, one can show that the effect of diffusion in 3 dimensions boils down to a perturbation of the form (14), with Q = P and E = k N , where k is the geometrical curvature of the scroll wave filament (i.e. the 3D extension of the rotation centres C) and N the local normal vector to it [19].
The resulting spiral and scroll wave dynamics can for both applications mentioned above be found by substituting (14) into the general law of motion (13). Here, we will assume that the RFs W M are essentially localized (as shown numerically in [25][26][27]) within an area of size d and that the spatial scale over which the fields E vary is larger than d. In the lab frame of reference, this delivers: In the case where E = k N , Q = P, system (15) describes the evolution of the scroll wave filament position X b in every plane locally orthogonal to the filament. We have thus generalized Keener's law of motion [19]. The main difference is that the coefficients Q M A depend on the meander phase ψ. The law of motion (15) can be simplified considerably by averaging it over several meander periods. This is most easily seen using Fourier series: where c.c. is the complex conjugate. The dynamics of the center of the meander flower in each cross-section perpendicular to the filament is then We note that unless in resonance, the set {| ω + kΩ| ∈ {0, ±1, ±2}, k ∈ Z} has a strictly positive minimal element, say ω min . Then, all non-constant terms in (17) will oscillate at a frequency of at least ω min .
However, the sole constant term F 0,0 b in the righthand-side of (17) will induce a constant drift velocity: In vector notation, this result can be written as where V is the net drift motion of the filament, T is the unit tangent to the filament for 3D scroll waves, and T = e z for 2D spiral waves in the XY-plane. From Eq. (18), the drift components parallel and perpendicular to the applied external field E are given by The time-averaged equation of motion for meandering spiral waves (19) exhibits the same dynamics as in the circular-core case. If h describes diffusive coupling in the third spatial dimension ( E = k N , Q = P), Eq. (19) happens to reduce to the circular-core result from [20]: where N and B are the normal and binormal vectors to the filament. Then, we can interpret Γ 1 and Γ 2 as the scalar and pseudoscalar filament tension. Since Eqs. (19) are the laws of motion for the filament of a meandering scroll wave, it follows from [20] that the period-averaged filament length increases monotonically in time if Γ 1 < 0 and decreases if Γ 1 > 0.
To validate our results, we have determined the coefficients P M A (ψ) for the Barkley and FK kinetics by applying E for a short time interval at different values of the meander phase ψ, see Supp. B for details of numerics. Averaging P B A (ψ) over one period delivered Γ 1 = −3.97, Γ 2 = 0.70 for Barkley kinetics and Γ 1 = 0.455, Γ 2 = 0.302 for FK kinetics. Theoretical predictions (19), (21) using the measured Γ 1,2 are in good agreement with the observed drift of spirals in constant field E, and with circular scroll ring dynamics, see Fig. 3.
Since the chosen parameters in Barkley kinetics yield Γ 1 < 0, the filament will undergo Euler buckling beyond a critical thickness, as we have already seen Fig. 2a. The FK model has Γ 1 > 0, and Fig. 2b shows that a transmural filament indeed relaxes to the minimal length.
Until now, it was assumed that perturbations are small and ω is not. As noted already in [39,40], if either condition is broken, phase-locking between spiral rotation and its meandering may happen. We are now in a position to describe this phenomenon quantitatively. For Barkley kinetics as in Fig. 1a where ω = 0.08 Ω = 1.25, one finds a qualitatively different tip trajectory when E = | E| ≥ 0.04, see Fig. 4a. From the first of Eqs. (15), one can show similarly to [41] that a necessary condition for locking the rotation phase is The locked rotation angle will be φ = arccos (−ω/Q) − arctan Q φ 2 /Q φ 1 . Given the computed Q φ A , expression (22) predicts E crit = 0.041, closely matching the value of 0.04 found in Fig. 4a. Fig. 4b shows a comparison for different values of parameter a; it can be seen that the Arnold tongue for phase-locking is well described by Eqs. (22).
In the (a, b) parameter space of Barkley's model, phase-locking is found near the line of resonant meander. Already for a field strength of E = 0.03, Fig. 4c shows phase-locking in a significant portion of the meander region, where it leads to relatively large drift velocities (see Fig. 4d). In qualitative terms, Fig. 4a shows that the meander flower opens up during phase-locking, and the resulting drift speed is therefore close to the mean 'orbital velocity' ωR of the tip along the meander flower, where R is the time-averaged radius of the meander flower. This result does not contradict Eq. (15) since when ω → 0, the centre of the rotating frame is far away. One can instead use a different rotating frame, with the origin shifted to the average tip position. This gives, in leading order, (23) We noted that the curve V = 0 closely matches the locus of resonant meander. We have however not found analytical proof of this property, and a counter-example in the Luo-Rudy-I cardiac tissue model is known [42].
Discussion. One motivation for this work was to see how the concept of filament tension generalizes to meandering scroll waves. Only when averaged over many meander periods, the dynamics reduces to the circular-core case. The tension concept has already been used for meander in cardiological literature [13] and modelling studies [42]. Here, we have shown that the emerging property of filament tension does indeed explain the (in)stability of scroll waves in simple cases, such as those in Fig. 2. Real heart tissue is more complicated in many respects. For instance, a significant phenomenon is pinning to heterogeneities. Some pinning effects have been described before using perturbation methods [43][44][45]. Thin domains of irregular thickness L(x, y) (e.g. the cardiac wall) can also be captured by (15), with h = P∇ ln L [45].
On short time-scales, dynamics is much more complex and the concept of filament tension cannot be applied. The orientation of the meander pattern may phase-lock to external fields, thickness or parameter gradients.
In general, the theory that was presented here opens the pathway to analysing and predicting the trajectory and stability of meandering spiral and scroll waves in reaction-diffusion media of diverse nature.
HD was funded by FWO-Flanders during part of this work. The computational resources (Stevin Supercomputer Infrastructure) and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by Ghent University, FWO and the Flemish Government department EWI. IVB and VNB gratefully acknowledge EPSRC (UK) support via grant EP/D074789/1. IVB acknowledges EPSRC (UK) support via grant EP/P008690/1. VNB acknowledges EPSRC (UK) current support via grant EP/N014391/1 (UK).
We integrate Eq. (1) forward in time by explicit Euler stepping with time step dt on a finite differences grid with spatial resolution dx and size N x × N y × N z . Diffusion is implemented using finite differences with a 5-point Laplacian in 2D and 7 points in 3D.
For Barkley kinetics, we used parameter values a = 0.58, b = 0.05, c = 0.02 unless stated otherwise. The model has dimensionless space and time units; we use N x = 500, dx = 0.1 and dt = 0.002375 in 2D and dt = 0.0016 in 3D. The spiral tip was found every 0.1 time units as the intersection of the isolines u = 0.5, v = 0.5a − b [30] using the algorithm in [17].
With the Fenton-Karma cardiac tissue model [17], we selected the guinea pig (GP) set of model parameters [17], since it yields a quasi-periodic meandering spiral with linear core. Notably, the reaction kinetic functions F(u) are not continuously differentiable with respect to u, so the Jacobian matrix F (u 0 ) in Eq. (4) contains singular contributions. However, in this study we do not directly solve the linearized equations, and F(u) can be regularised by approximating Heaviside function in its definition with a smooth sigmoidal function, with any required accuracy. We do not carry out this limit procedure here, but the predictions of our theory agree well with the outcome of numerical simulations from the unmodified Fenton-Karma GP model. We used P = diag(0.1, 0, 0) mm 2 /ms as in [17] and we decreased the spatial grid to dx = 0.15 mm to find quasi-stationary rotation. A time step of dt = 0.053 ms was chosen in a grid of size 400 × 400. The tip line was tracked as the intersection of the isosurfaces u = 0.5 and ∂ t u = 0, resulting in the tip trajectory of Fig. 1b  If the response functions W M are known with sufficient accuracy, the matrix elements Q M N can be found by evaluating the integrals, in a manner similar to the circular-core case [36? ]. In practice, however, it is simpler to measure the drift induced by an external field which is imposed during a given fraction of the meander cycle. First, we measure the absolute position and phase of an unperturbed meandering spiral wave as detailed in [27] and denote it as X M 0 = (X 0 , Y 0 , φ 0 , ψ 0 ). Thereafter, we run many simulations that deliver a global stimulus at different phases of the meander cycle: for different t p we apply a time-dependent field E(t) = E 0 H(|t − t p | − ∆ t /2) e x of duration ∆ t , where H is the Heaviside step function. For each simulation, we measure the absolute spiral phases and position as explained in [27] in a time interval [t start , t end ], where t start t p . The length of the interval was taken to be 120 ≈ 24 T for Barkley kinetics and t 2 = 1 s ≈ 40 T for the FK guinea pig model. By comparing with a reference case without perturbation, we can compute the matrix elements: We repeat the procedure for a field along the y-direction to find Q m y (t p ) and then convert the matrix elements to the periodic functions Q m A = Q m a R a A (φ) which only depend on ψ. Note that 2π-periodicity of these functions is not guaranteed by construction, but is indeed observed with good accuracy in Fig. 5, which we take as an indication that the results are reliable and the analysed spiral wave regimes are indeed biperiodic. Fig. 5 shows the computed curves Q m A for the Barkley and FK models with parameters as above. Since we took Q = P, the computed coefficients are denoted P m A . For the cardiac tissue FK-model, these results should be interpreted in terms of filament tension rather than electrophoretic drift, which is reserved for chemical systems. To find a suitable E 0 we tried the method first at ψ = 0 for various field strengths and chose the maximal E in the model for which the response was still linear in the field strength. This resulted in E = 0.5 for Barkley kinetics and and E = 0.1 mm −1 for the FK model.