A response function framework for the dynamics of meandering or large-core spiral waves and modulated traveling waves
Hans Dierckx, Alexander V. Panfilov, Henri Verschelde, Vadim N. Biktashev, Irina V. Biktasheva
AA response function framework for the dynamics ofmeandering or large-core spiral waves and modulated traveling waves
H. Dierckx, A. V. Panfilov,
1, 2
H. Verschelde, V.N. Biktashev, and I.V. Biktasheva
4, 3 Department of Physics and Astronomy, Ghent University, 9000 Ghent, Belgium Laboratory of Computational Biology and Medicine, Ural Federal University, Ekaterinburg, Russia College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4 4QF, UK Department of Computer Science, University of Liverpool, Liverpool L69 3BX, UK (Dated: January 25, 2019)In many oscillatory or excitable systems, dynamical patterns emerge which are stationary orperiodic in a moving frame of reference. Examples include traveling waves or spiral waves in chemicalsystems or cardiac tissue. We present a unified theoretical framework for the drift of such patternsunder small external perturbations, in terms of overlap integrals between the perturbation andthe adjoint critical eigenfunctions of the linearised operator (i.e. ‘response functions’). For spiralwaves, the finite radius of the spiral tip trajectory as well as spiral wave meander are taken intoaccount. Different coordinates systems can be chosen, depending on whether one wants to predictthe motion of the spiral wave tip, the time-averaged tip path, or the center of the meander flower.The framework is applied to analyse the drift of a meandering spiral wave in a constant externalfield in different regimes.
I. INTRODUCTION
Spiral waves are remarkable self-sustained patternswhich arise in various extended systems, including oscil-lating chemical reactions [1], catalytic oxidation [2] andbiological systems. In biology, spiral waves have beenobserved across vastly different spatial scales, as they or-ganise intracellular calcium waves [3], slime mould ag-gregation [4], waves of spreading depression in the brain[5], and cardiac contraction during arrhythmic events[6–9]. In the context of cardiac arrhythmias, differentdrift regimes are thought to correspond to different heartrhythm disorders [7]. Moreover, tracking of rotation cen-ters of cardiac rotors was shown to improve ablation ther-apies of atrial arrhythmias [10].The development of asymptotic description of drift ofspiral waves has been contingent on two important issues:localization of adjoint symmetry modes [11–16] and anappropirate choice of the drift system of reference [17–20]. On this pathway, the choice of an “optimal” systemof reference naturally depends on the unpertubed spiralwave dynamics, i.e. whether it is (i) a rigidly rotatingspiral wave with its tip following a perfect circle, (ii) abi-periodic “meander” regime when the spiral wave pe-riodically changes its shape and its tip describes a bi-periodic flower-like trajectory, or (iii) more complicatedcases dubbed ‘hypermeander’. While the asymptotic the-ory of the drift of rigidly rotating spiral waves is the mostadvanced [19, 20], the theory of hypermeander drift is theone least developed [21–23]. The current work aims tolay out a technical framework for the drift of a bi-periodicmeander regime, and illustrate its application on a simpleexample of electroforetic drift.The propagation of cardiac excitation, as well as manyother systems in which spiral waves are observed, canbe described by reaction-diffusion (RD) systems (see, e.g. [24]), ∂ t u = ∂ j [ P jk (⃗ r, t ) ∂ k u ] + F ( u , ⃗ r, t ) . Here, u = u (⃗ r, t ) ∈ R n is the list (column-vector) of n state variables, u = ( u , . . . , u n ) T , e.g. the concentra-tions of reacting species, ⃗ r = ( x j ) ∈ D ⊂ R N is thevector of Cartesian coordinates, N = , t is time. The vector-function F ∶ R n × R N + → R n de-scribes the local dynamics of state variables, e.g. the re-action rates. The space-time dependent matrix-tensor ( P jk ) ∶ R N + → R n × n × R N × N contains the diffusivitiesof the state variables on its diagonal (of which some maybe vanishing). Non-diagonal terms of P are present insystems with cross-diffusion; see e.g. Ref. [25] and Fig.1(b).The asymptotic description follows the paradigm thatfirst an idealized mathematical situation is considered,which is an unbounded, perfectly stationary, homoge-neous and isotropic medium governed by the isotropicreaction-diffusion equation (RDE) ∂ t u = P ∆ u + F ( u ) (1)and then any deviations from the idealized picture areconsidered perturbatively. Here the diffusion matrix P is constant and the kinetic functions F ( u ) are taken uni-form in space and time. In the context of pattern forma-tion, a bistable, oscillatory or excitable point system iscommonly used. For well-chosen sets of reaction kinetics,i.e. functions F ( u ) , the system (1) may sustain waves offinite amplitude, which propagate through the medium,see Fig. 1(a–b).Although our methodology works for any reaction-diffusion system allowing pattern formation, we will usebelow examples with the reaction kinetics of FitzHugh- a r X i v : . [ n li n . PS ] J a n Nagumo [14, 26] P = diag ( , ) , F ( u ) = [ α − ( u − u / − u ) , α ( u − γu + β )] T , Barkley [27] P = diag ( , ) , F ( u ) = [ ς − u ( − u )[ u − ( u + b )/ a ] , u − u ] T and Fenton and Karma. The equations and parametersof the latter model (with 3 state variables) are found in[28].In N ≥ u j ( x, y, t ) = u ∗ j , u k ( x, y, t ) = u ∗ k . (2)where j , k are indices of two selected state variables, j ≠ k , and u ∗ j , u ∗ k are two appropriately chosen constants.In N = j = , u ∗ j = . , k = , u ∗ = . a − b . It is often alsoconvenient to take ∂ t u j = χ in a suitably chosen plane in the phase space of localkinetics: tan χ (⃗ r, t ) = u k (⃗ r, t ) − u ∗ k u j (⃗ r, t ) − u ∗ j , (3)This mapping will produce a phase singularity (PS) inthe vicinity of the spiral tip [30]. If j, k, u ∗ j , u ∗ k are setto the same values in Eqs. (2) and (3), the PS coincideswith the spiral wave tip definition. Otherwise, differentobservers will generally not exactly agree on the tip orphase singularity position [31].In three spatial dimensions, the solution to thereaction-diffusion equation (1) consisting of a stack of spi-ral waves is known as a ‘scroll wave’. When tracking thescroll wave’s tip or PS, Eqs. (2) and (3) produce a curvethat is known as the ‘scroll wave filament’. The same re-mark as in 2D holds here: different algorithms and chosenthresholds will yield different filament curves, which gen-erally lie in each other’s vicinity, i.e. in the tubular scrollwave core region.However, if the tip trajectory is circular, all tip pathsand PS trajectories will in 2D describe circles around a unique rotation center C at the same angular frequency ω . In previous works, the dynamics of circular-core spi-ral waves have been analysed in terms of the motion ofthe rotation center C [32–35]. For circular tip paths, theregion inside it never excites and is sometimes called the‘spiral wave core’. The dynamics of circular-core spiralwaves is reasonably well understood, and their drift re-sponse to small external stimuli can be found by project-ing the stimulus on the inertial manifold of the systemusing so-called ‘response functions’ [17, 19, 36, 37].Interestingly, in several experiments and numericalsimulations, both in chemical [38, 39] and biological[9, 28, 40–44] systems, the spiral wave was found to per-form not a rigid rotation, but a more complicated motion,a phenomenon known as ‘spiral wave meander’ [45–47].In excitable models, meander can arise due to the wavefront interacting with the wave back from the previousexcitation. In the simpler cases, the meander is quasi-periodic, in the sense that the evolution is periodic upto a Euclidean transformation. This transformation canalways be written as a pure rotation or a pure translation(see below), which is why one can refer to it as ‘biperi-odic’ meander, even when the tip trajectory is not a su-perposition of two circular motions. In more complicatedcases, often called “hypermeander” [47], the tip motioncan have more than two periods [22] or be chaotic, result-ing in deterministic Brownian motion of the tip position[21]; we will not consider such cases here.A significant step in understanding meander of spiralwaves was made by Barkley et al. [48] and Karma [49],who showed that the transition to meander from thecircular-core case exhibits typical features of a Hopf bi-furcation. Barkley [50] then showed that the process canbe described by a set of five coupled ordinary differentialequations. As a result, a tip trajectory is formed that is asuperposition of two circular rotations, making a flower-like tip trajectory. This picture was put on constructivefooting in [51–54], using “skew-product” representationof the reaction-diffusion system exploiting its symmetrywith respect to Euclidean motions of the plane. Theseworks described the equivariant Hopf bifurcation (i.e. theregime close to the transition to meander), while manychemical and cardiac models show a qualitatively differ-ent tip trajectory, i.e. zig-zagged, star-like path knownalso as the “linear core” case. It turned out, however,that skew-product representation does not need to be re-stricted to the vicinity of the transition zone, and themore complicated meander pattern can be described asrelative periodic solutions not necessarily related to aHopf bifurcation [55, 56].Correspondingly, evolution of biperiodically meander-ing spirals, both flower-like and star-like, in response tosmall perturbations, can be analysed by linearising thesystem around a relative periodic orbit [16]. To showthe steps of this process in detail is the purpose of thismanuscript; it extends the procedure for rigidly rotatingspirals found in e.g. [14, 19, 36]. Since periodically de-forming waves in one spatial dimension are obtained as a) b) c)d) e) f)FIG. 1. Qualitatively different one- and two-dimensional patterns that arise as solutions to the reaction-diffusion equation (1),depending on the chosel reaction model F ( u ) and initial conditions. a) Stationary traveling wave in the FitzHugh-Nagumomodel [26] ( α = . β = . γ = . t = t = a = . b = . ς = . a = . b = . ς = . a = . b = . ς = . a special case, we also describe their dynamics.A significant part of this paper will be dealing with in-troducing different coordinate systems that are suitableto capture the drift dynamics of solutions to the RDE(1) because, depending on which definition of filamentone adopts one can obtain simple or complex laws ofmotion. These laws can be further simplified if one av-erages in time over rotation cycles of the spiral or scrollwaves [36]. For non-stationary spiral or scroll waves, apossible averaging method is to analyse the motion ina Fourier series and only keep the non-oscillating, sec-ularly growing terms [57, 58]. Some recent works havecomputed the leading order eigenmodes in the labora-tory frame of reference [15, 59]. At the end, differentdescriptions need to be compatible, and it is our aim toillustrate that they all describe the same dynamics.This paper is organised as follows. In Sec. II we dis-cuss the concept of symmetry breaking by the formedpatterns, which will determine the number of degrees offreedom required to uniquely describe such patterns. Wefirst illustrate this concept in Sec. III on traveling wavesin one dimension, treating the cases of rigidly movingand periodically deforming waves. In Sec. IV we con-sider spiral wave patterns in two spatial dimensions, inthe regimes of rigid rotation, meander and resonant me-ander. We consider in detail multiple possible frames of reference, which are defined with respect to either instan-taneous or time-averaged dynamics. In Sec. V we lin-earise around the solution and discuss the critical eigen-modes of the system. This knowledge is used in Sec. VIto derive the equations of motion for modulated wavepatterns, in a manner that is valid for both traveling andspiral waves. In Sec. VII we re-interpret our approach interms of the geometry of the phase space of the reaction-diffusion system considerd as a dynamical system. Thispart is optional, but we found it useful to include as fewworks make the connection between the ‘physics’ styleof description and the more abstract dynamical systemsviewpoint.In Sec. VIII the drift response of a spiral wave in aconstant external vector field is analysed. We show howthe dynamics can be averaged over rotation and temporalphases in order to find a simpler, manageable equationof motion. In particular, we demonstrate that if a mean-dering spiral wave is phase-locked to a constant externalfield of magnitude E , its drift velocity can be larger thanorder E : it is close to the ‘orbital velocity’ which whichthe spiral wave circumscribes the meander flower.We conclude with a discussion of our present work inSec. IX. To assist the reader in remembering the nota-tions used, we provide a list of notations in the Appendix. II. SYMMETRY CONSIDERATIONS
It is customary to consider Eq. (1) in a bounded spa-tial domain with Neumann boundary conditions, whichagrees with many physical situation where such modelsare relevant, and this is used in all numerical simula-tions of this work. However, to develop the theory ofsymmetry-breaking, we here consider the whole space R N , with N = N = W M defined in Eqs. (17) belowdecay exponentially far from the wave front (in 1 spatialdimension) or far from the spiral wave tip (in 2 spatialdimensions). Therefore, faraway boundaries are not feltat e.g. the wave front or spiral tip.Since the system (1) involves the spatial variables onlyvia the Laplacian ∆ in the whole space R N , it is equiv-ariant (symmetric) with respect to the isometric trans-formations γ ∶ R N → R N of this space, meaning that if u (⃗ r, t ) is a solution, then ˜ u (⃗ r, t ) = u ( γ − ⃗ r, t ) is also a solu-tion, for any isometry γ . For our purposes it is sufficientto consider only the orientation-preserving transforma-tions, that is, γ ∈ SE ( N ) . Besides, being autonomous,this system is equivariant with respect to the group oftranslations of the time axis, SE ( ) ∼ R , that is if u (⃗ r, t ) is a solution, then ˜ u (⃗ r, t ) = u (⃗ r, t − T ) is also a solutionfor any T ∈ R . (Note that this is different from sayingthe solution is ‘invariant’ with respect to time transla-tion, which would imply a solution that is constant intime.) Thus, we can say overall that (1) is symmet-ric with respect to the inhomogeneous Galilean group G = SE ( N ) × R . The dimensionality of this Lie group isdim ( G ) = N ( N + )/ + u (⃗ r, t ) of Eq. (1) may be invariantwith respect to none, or all, or part of the symmetriesof the equation itself. Invariance of the solution withrespect to a purely spatial symmetry γ ∈ SE ( N ) can ingeneral be reduced to a lower-dimensional solution and isnot further considered. E.g. a plane wave in N =
2, whichis invariant for translations along the wave front. Thiscase can be reduced to the study of a traveling wave in N = N = c , u (⃗ r + cT , t + T ) = u (⃗ r, t ) forany T , which implies that there are members of G thatleave the solution invariant.Let the group of symmetries leaving a solution invari-ant be H , H ≤ G . If this is also a Lie group, then ap-plication of all possible transformations γ ∈ G to u (⃗ r, t ) will produce a manifold of solutions, of dimensionality N BS = dim G − dim H . We refer to this dimensionalityas the number of broken symmetries (implying that the “unbroken ones” are in the subgroup H ). By construc-tion, the tangent vectors to this manifold of solutions aresolutions of the (1) linearised on u (⃗ r, t ) which do not ex-ponentially grow or decay with time, so they make thecenter subspace of this solution. If the solution is an at-tractor to the system, this manifold of solutions will becalled an ‘inertial manifold’, see details in Sec. VII.In the following text, it will sometimes be convenientto indicate for every quantity, e.g. wave front or tip coor-dinates, the reference frame in which it is defined, usingdifferent math accents. We shall call the ‘simplest’ framethat makes the solution constant or periodic the centerframe of reference . For any quantity f , its value pertain-ing to this frame of reference, will be denoted as ○ f . Theframe which best follows the wave front or spiral wavetip will be called the ’tip frame‘, with the correspondingnotation ˇ f . An list of all notations and symbols usedthroughout this manuscript is given in the appendix. III. ONE-DIMENSIONAL TRAVELING WAVESA. General framework
Let us first consider a generic choice for the frame ofreference to describe a traveling wave U ( x, t ) . We de-note the origin of the lab frame with Cartesian coordi-nate x as O . In the lab frame, we can define the wavefront position x = Z( t ) by thresholding a state variable u j (Z( t ) , t ) = u ∗ j and demanding that ∂ t u j (Z( t ) , t ) > Z( t ) , with lab frame coordinate x = Z ( t ) . However, al-ternative definitions are possible, e.g. for the wave shownin Fig. 1(b) the barycentre of ∣ u ∣ was used to denotewave position.Next, we introduce a (yet unspecified) moving frameof reference with origin X that has lab frame coordinate X ( t ) and denote the new spatial coordinate in the mov-ing frame as ρ : ρ = x − X ( t ) , t = τ . (4)The moving frame can be chosen in different manners;an overview is given in Tab. I.The instantaneous wave front position is denoted Z ( t ) ,it is found at the position ρ = ζ ( t ) in the moving frame: Z ( t ) = X ( t ) + ζ ( t ) . B. Symmetry breaking
In a one-dimensional homogenous spatial domain ( N = G ∼ R × R , and dim ( G ) =
2. A traveling wave type N BSS N BTS N BS Frame name Origin Choice Frame velocityany 1 0 or 1 1 or 2 lab fixed
X =
O c = X = Z c = c f periodically modulated front 1 1 2 fully co-moving at front X = Z c ( Ψ ) = c f ( Ψ ) constant-speed co-moving frame near front X ≠ Z c = c f ( Ψ ) = constTABLE I. Overview of reference frames used in the text to describe the motion of 1D traveling wave patterns. N BSS : numberof broken spatial continuous symmetries, N BTS : number of broken temporal continuous symmetries, N BS = N BSS + N BTS :number of broken continuous symmetries, O : lab frame origin, Z : instantaneous wavefront position, X : chosen origin of themoving frame, c f : speed of the wave front in the lab frame, Ψ: temporal phase describing the progression of the solution alonga deformation cycle (relative periodic orbit). traveling-wave solution does not have translational sym-metry, neither in time nor space, so we can say it breaksthe translational symmetry. However, if the wave travelswith constant speed c f without changing its shape, thenwe can say that u ( x − X, t − T ) = u ( x, t ) as long as X and T are related by X = c f T . That means, that the solu-tion is invariant with respect to a one-parametric groupof transformations, so dim H =
1, and N BS =
1. This isthe case of rigidly translating waves, to be treated in sub-section III C. However, if there is no frame of referencemoving with constant speed in which the solution is con-stant in time, then the subgroup H is trivial or at mostdiscrete, so dim H =
0, and N BS =
2. The case wherethis dynamics is periodic, so H consists of translations x → x + (cid:96)c f T , t → (cid:96)T , (cid:96) ∈ Z , is treated in subsection III D.The non-periodic case, when H is trivial, falls outsideour present scope. C. Rigidly translating waves
Firstly, if the wave is rigidly moving, Z ( t ) = Z ( ) + c f t with c f the constant traveling wave speed. In this case,we simply take X ≡ Z and X ≡ Z and speak of the ‘co-moving frame’. In this frame, the traveling wave is asolution to the RDE (1) of the form: U ( x, t ) = u ∗ ( x − X ( t )) , (5a) ∂ τ X = c. (5b)Note that Eq. (5a) describes the type of solution and(5b) gives the equation of motion for the collective coor-dinate X ( t ) (wave front position) that is introduced sincethe pattern breaks the translational symmetry along x .If a perturbation of order η is added to the system, Eqs.(5a) and (5b) will gain additional terms of O( η ) in theirright-hand sides, describing the wave profile change andwave speed correction, respectively. D. Periodically modulated waves
Traveling wave may also possess a variable shape andvariable propagation speed c f , in which case X ( t ) = X ( ) + ∫ t c f ( t ′ ) d t ′ . This situation has been encountered both in experiment and simulations [25, 60], see Fig. 1bfor an example. We only treat periodically modulatedwaves here, where c f ( t ) has temporal period T . In thepresence of perturbations, it will be convenient to rescalethe period to 2 π , define the positive constant Ω = π / T and call Ψ = Ω t the temporal phase of the solution, whichlabels how far the wave front shape has gone through itscycle. That is, Ψ is essentially a scaled time variable, andΩ is the scaling constant.Throughout this work, we will use the bar notation todenote the average over one temporal cycle: f = π ∫ π f ( Ψ ) dΨfor any 2 π -periodic f ( Ψ ) . We now have multiple optionsto choose the frame of reference.The traveling wave solution is now periodic in a movingframe, i.e. U ( x, t ) = u ∗ ( x − X ( t ) , Ψ ( t )) ,∂ t X = c ( Ψ ) ,∂ t Ψ = Ω . (6)Due to the scaling of time, u ∗ is 2 π -periodic in Ψ.
1. Fully co-moving frame
In the ‘front frame’, we take again
X ≡ Z and X ≡ Z , atthe expense of a non-constant propagation speed c f ( Ψ ) .
2. Constant-speed co-moving frame
Alternatively, one may choose to let the frame move atthe time-averaged velocity c = c f ( Ψ ) . Since c is constant here, this frame is simpler than thefully co-moving frame. Keep in mind, however, that X now describes the average progression of the wave frontposition; the true wave front position Z can be foundfrom Z = X + ζ, ζ ( Ψ ) = ∫ Ψ0 ( c f ( Ψ ′ ) − c ) dΨ ′ . IV. SPIRAL WAVESA. General framework
Let x a , a ∈ { , } be Cartesian coordinates in the labframe with origin O . In the lab frame, the spiral wavetip Z as defined by (2) describes a path x a = Z a ( t ) thatmay be circular or flower-like, see Fig. 1.To capture spiral wave dynamics, we will introduce amoving frame with origin X ( t ) , whose lab frame coor-dinates are x a = X a ( t ) . At all times, the orientation ofthe new frame of reference with respect to the lab frameis given by the rotation phase Φ ( t ) , with rotation rate ∂ t Φ = ω that may be time-dependent. So, ω takes posi-tive values for counterclockwise frame rotation and nega-tive values for clockwise frame rotation. If ω is constant,we introduce K = sgn ( ω ) ; if it is periodic, K = sgn ( ω ) .Remark that for meandering spirals, ω coincides with thespiral rotation rate only if the ‘tip frame’ (defined below)is chosen. For the other frame choices in Tab. II, ω equalsthe (typically slow) precession rate of the meander pat-tern.The spatial coordinates in the moving frame are de-noted ρ A , and co-moving time coordinate is τ . Thegeneric coordinate transformation between the frames isthus given by: x a = X a + R aA ( Φ ) ρ A , t = τ . The rotation matrices connecting lab frame coordinates x a to rotating coordinates ρ A are given by ( R Aa ( Φ ) ) = ( cos Φ sin Φ − sin Φ cos Φ ) , ( R aA ( Φ ) ) = ( cos Φ − sin Φsin Φ cos Φ ) . Since the generator of rotations is the Levi-Civita sym-bol (cid:15) AB ( (cid:15) = − (cid:15) = ∂ Φ R Aa ( Φ ) = − R Aa ( Φ ) (cid:15) ab = − (cid:15) AB R Ba ( Φ ) ,∂ Φ R aA ( Φ ) = − R aA ( Φ ) (cid:15) AB = − (cid:15) ab R bA ( Φ ) . We will write the polar angle in the chosen frame of ref-erence as θ , while the spatial orientation of the referenceframe will be the rotation phase Φ. Due to rotationalinvariance of the RDE (1), the solution will generally de-pend only on ϑ = θ − Φ . This relation is the angular equivalent of Eq. (4).In the moving frame, the spiral wave tip position isfound at position ρ A = ζ A ( τ ) . Therefore, the lab frametip position Z a is given by: Z a = X a + R aA ( Φ ) ζ A . Differentiating with respect to time delivers ∂ t Z a = ∂ t X a + R aA ( Φ ) Ω ∂ Ψ ζ A − ω(cid:15) ab R bA ( Φ ) ζ A . a) b)FIG. 2. Different coordinate systems suitable to analyse andcircular-core spiral wave. a) Center frame: the origin is takenat the center of the circular tip trajectory. b) Tip frame: theorigin is taken at the current position of the spiral wave tip,and the ρ axis is tangent to the tip trajectory. B. Symmetry breaking
The Euclidean plane R is invariant under the Eu-clidean group SE ( ) of translations in x and y and rota-tions. However, the spiral wave solutions are not invari-ant under SE ( ) , since shifting or rotating a solution willyield a solution to the RDE (1) that is distinct from theoriginal one. So, at least 3 symmetries are broken by thespiral wave solution.As in the case of 1D pulses, one may introduce a mov-ing frame of reference, which in this case is also rotating;see below for different options. If in that frame, the spiralwave is stationary, one has N BS = N BS =
4. Here, we consider onlythe case where the spiral wave solution in the co-movingframe of reference is periodic in time. Returning to thelab frame of reference, this case represents a meanderingspiral wave, either in the non-resonant (see Par. IV D)or resonant case (see Par. IV E).
C. Rigidly rotating spiral waves
When going to a frame of reference rotating at ω = Kω s , (i.e. the spiral wave frequency), the solution be-comes periodic: U ( x a , t ) = u ∗ ( ρ A ) ,∂ t X A = c A ,∂ t Φ = ω (7)with c A and ω constant. We distinguish two frames ofreference here:
1. Center frame
The simplest coordinate system is found when takingits origin in the middle of the circular tip trajectory, as
Spiral type N BSS N BTS
Frame name Origin of moving frame Choice Rotation Translationany 3 0 or 1 lab fixed
X =
O ω = c a = X = C ω = Kω s c a = X = Z ω = Kω s c A = const ≠ X = C ω = β / T c a = X ≠ C , Z ω = β / T c A = const ≠ X ≠ C , Z ω = β s / T c A = const ≠ X = Z ω ( Ψ ) = β s / T c A ( Ψ ) resonantmeander 3 1 minimally rotating co-moving frame on a line X ≠ Z ω = c A = const ≠ X ≠ Z ω = K Ω c a = const ≠ X = Z ω ( Ψ ) = K Ω c A ( Ψ ) = N BSS : numberof broken continuous spatial symmetries, N BTS : number of broken continuous temporal symmetries, O : lab frame origin, C :middle of the tip trajectory, Z : spiral wave tip, X : chosen origin of the moving frame (“filament”), c a : velocity of X in the labframe, c A : its velocity in the rotating frame, K = ± ω , introduced sothat ω s > >
0, respectively. N BSS + N BTS = N BS , the total number of broken symmetries. shown in Fig. 2(a), so that c A = . This frame was used in the vast majority of previousworks on spiral wave theory [17, 18, 32–35, 57, 61, 62].This approach, however, is not optimal for describingthe interaction of the spiral wave with localised hetero-geneities, as dynamics is most influenced by stimuli closeto the spiral wave tip rather than rotation centre, andthe distance between the two can be significant.
2. Tip frame
Therefore, we introduce another frame of reference,with the origin of coordinates ρ A at the spiral wave tip, X = Z as shown in Fig. 2b: c A ≠ . In this ‘tip frame’, the velocity of the spiral wave tip c A becomes constant. One can even choose the orientationof the frame such that c = c = ωr with r the radiusof the tip trajectory.The tip frame closely relates to the kinematic approachof Zykov et al. [63] and is particularly suited to describelarge-core spiral waves.Finally, note that the freedom in the choice of the def-inition of the tip may be exploited by taking u ∗ j = u j (C) , u ∗ k = u k (C) , so that Z = C ; then we have c A = D. Meandering spiral waves
The case of biperiodic meander studied in this workis characterized by the fact that the spiral wave solution can be made time-periodic in a moving frame of refer-ence. As in the case of periodically modulated waves, weintroduce a temporal phase Ψ, which in the absence ofperturbations increases from 0 to 2 π over the time inter-val T : Ψ ( t ) = Ω t , Ω = π / T >
0. Without loss of general-ity, we may take Ψ ( ) =
0. In the context of spiral waves,we refer to Ψ as the meander phase, or temporal phase,as opposed to the rotation phase (angle) Φ.A meandering spiral solution is thus of the form: U ( x a , t ) = u ∗ ( ρ A , Ψ ( t )) , (8a) ∂ t X A = c A ( Ψ ) , (8b) ∂ t Φ = ω ( Ψ ) , (8c) ∂ t Ψ = Ω . (8d)with u ∗ π -periodic in Ψ.An overview of reference frames discussed below isgiven in Tab. II.The Euclidean transformation mapping u ( x, y, t ) to u ( x, y, t + T ) can, without loss of generality, be alwaysthought of as either a pure rotation or a pure translation,since the composition of a rotation and a translation canalways be written as a single rotation in 2D. Below wefirst treat the ‘pure rotation’ case, corresponding to non-resonant meander; the second case (pure translation),corresponding to resonant meander will be discussed inSec. IV E.
1. Definition of the angle β between petals In the case of non-resonant meander, the spiral patternrepeats itself after time T , but rotated around an angle β around a point C that will become the center of thespiral tip trajectory (‘meander flower’); see Fig. 3.Since we only impose that the solution be periodic, β can be chosen up to an integer multiple of 2 π . Twovalues can play a special role, however. First, we can a) b) c)d) e) f)FIG. 3. Different coordinate systems suitable to analyse and predict the motion of a meandering spiral wave. Panels a-d showcycloidal meander in Barkley’s model with parameters as in Fig. 1d; panels e-f illustrate linear-core meander in the FK-GuineaPig cardiac model. Panels (a,e) show the ‘center frame’, where the origin is taken at the center of meander flower, with theframe rotating at the constant rate ω = β / T , where β is the angle between consecutively visited petals. Panels (b,c) illustrateframes tailored describe the finite-core, which may either minimally rotate (over an angle β (see b) or co-rotate with the spiral,i.e. rotate over β s = β + πK where K = ± X of thecoordinate system at the current tip position and the ρ -axis tangent to the tip trajectory. The rotation rate of the tip frameis ω ( Ψ ) , non-constant but periodic. define this angle by following the movement of the spiral.If we attach a reference direction to the spiral wave tip(e.g. ⃗∇ u evaluated at the tip), this vector will turn overthe angle β s during one spiral period, and for a farawayobserver, the average spiral wave frequency (expressed asa positive number) will be ω s = Kβ s / T , where K = ± K = K = − β as the element of theset β s + π Z that has the minimal absolute value, andcall this minimal value β . For the cases studied here,we find β s = β + πK with 0 < ∣ β ∣ < π .
2. Center frame
We pick a frame rotating at constant frequency ω = β / T around C , i.e. we take X = C . In this ‘center frame of reference’, the meandering spiral solution becomes pe-riodic. The solution is given by (8), with c A = , ω = β / T .
This set of equations describes the spatial position thespiral wave in terms of the position of the center of themeander flower, and was used in [16]. Since we chose
X = C , Eq. (8b) expresses that the center of the meanderflower does not move in the absence of perturbations tothe RDE (1).It is also possible to relax the convention to take the pe-riod T as the minimal time interval in which the solutionis periodic modulo Euclidean group actions. E.g. whenlinear cores rotate over β ≈ ○ [28], one could define theperiod to span two such cycles, bringing β ≈ ○ instead.This formalism can then be used to study phase-lockingto constant external fields, generalising the results in [16].Another theoretical possibility is to use the same set-ting as described by system (8), but with the choice β = β s instead of β . This choice would ensure that thespiral wave solution, considered in the rotating frame, os-cillates but not rotates and remains approximately sta-tionary far from the centre. Neither of these two theoret-ical alternatives is used later in this paper, so we mentionthem only for completeness.
3. Minimally rotating finite-core frame
It can be advantageous to allow the origin of the coor-dinate system to better follow the spiral wave tip. Onepossible choice is shown in Fig. 3b for reaction kineticswith cycloidal meander. We let the new origin move atconstant angular velocity on a circular trajectory thatapproximates the exact tip trajectory, and refer to thissituation as the ‘minimally rotating finite-core frame’: ∂ τ c A = , ω = β / T .
4. Co-rotating finite-core frame
We can let the minimally rotating finite core framerotate around its axis, at a rate K Ω, to better follow thespiral wave rotation, bringing the absolute rotation rateto β s ; see Fig. 3c. We will refer to this case as the ‘co-rotating finite-core frame’. It is again described by Eqs.(8), with ∂ τ c A = , ω = β s / T . (9)In this frame, the spiral solution u ∗ ( ρ A , Ψ ) will notrotate if Ψ is increased from 0 to 2 π .
5. Tip frame
Lastly, we let the origin X of the co-moving frame ofreference coincide with the spiral wave tip position, i.e. X = Z , X a ( t ) ≡ Z a ( t ) ; see Fig. 3c and f. Then c A and ω become non-constant, 2 π -periodic functions of Ψ.After one meander period of duration T , the frame willhave turned over an angle β s = ∫ T ω ( Ψ ( t )) d t = ∫ π ω ( Ψ ) dΨ = ¯ ωT . If the tip trajectory is smooth, we can without lossof generality suppose that c = c = c ( Ψ ) . This canalways be achieved by exploiting the freedom in the def-inition of frame orientation angle Φ.
6. Comparison with the classical theory of meander
In the classical theory of rigidly rotating spiral waves,the critical eigenmodes for translation have eigenvalues ± iω , where ω is the rotation frequency of the rigidlyrotating spiral wave. If due to a parameter change inthe medium, another complex eigenvalue pair crosses theimaginary axis at ± iω , we have a Hopf bifurcation inthe quotient system (rotating frame). If this bifurcation is supercritical, this results in an epicycloidal tip trajec-tory as shown in Fig. 3(a–c). In Barkley’s notation [50],it is seen that beyond the Hopf bifurcation the absolutetip velocity will oscillate at ω . Hence, we conclude thatBarkley’s ω equals Ω used throughout this work. Sec-ondly, in the circular-core regime just before the Hopfbifurcation, the spiral’s rotation frequency is ω , corre-sponding to ω s in our notation. We conclude that in thetip frame (where β = β s ), ω = ω and in the minimallyrotating frames, (where β = β ), ω = ω − K Ω.The case of resonant meander happens in the classicaltheory when ω → ω . In our description, this happenswhen ω → ω → Ω in the tip frame or fully co-rotating frame.
E. Resonant meander
We now consider the case where the meandering spiralreturns after the temporal period T to the same state, buttranslated over a vector d a = c a T ; see Fig. 1(f). Differentauthors call this sort of solutions cycloidal motion [46], 0 ○ isogon contours [47], modulated travelling waves [50, 64],resonant linear motion [53] or linear drift [23]. We shallrefer to it here as resonant meander.Note that during dynamics under external perturba-tions, the direction of spontaneous drift may change, suchthat a rotation phase needs to be introduced neverthe-less. This makes this case distinct from the 1D modu-lated traveling wave from Par. III D.The solution is still given by (8), with c A and ω de-pending on the chosen reference frame:
1. Minimally rotating co-moving frame
First, we let the origin of our coordinate system moveon a straight line with velocity c , c , without frame rota-tion, parallel to the line along which the resonant spiraltravels over time. This situation is sketched in Fig. 4a.Then, we find: ∂ τ c A = , ∂ τ Φ = . Since the frame is not rotating, the function u ∗ ( ρ A , Ψ ) shows a spiral that rotates a full turn when Ψ is raisedfrom 0 to 2 π .
2. Co-rotating, co-moving frame
We can also let the previous frame rotate at ω = K Ω,while its origin moves on a straight line at constant speedparallel to the displacement direction of the resonant spi-ral. We again obtain Eqs. (8), with c A = R Aa ( Φ ) c a , ∂ τ Φ = K Ω . a) b)c)FIG. 4. Reference frames for resonant meander. a) The ‘min-imally rotating co-moving frame’ is rigidly translating suchthat the solution becomes periodic. During one period (i.e.between the frames shown), the frame rotates over an angle2 πK . b) The ‘co-rotating co-moving frame’ is obtained fromthe average frame by letting the moving coordinates rotateat constant frequency K Ω = πK / T . c) Tip-based frame ofreference, with the origin of the coordinate system at the cur-rent tip position and the ρ -axis tangent to the tip trajectory.Rotation rate of the frame is ω ( t ) , non-constant but periodic. and c a constant. The spiral wave profile u ∗ ( ρ A , , Ψ ) does not rotate a full turn when Ψ is raised from 0 to 2 π ;it only deforms (‘breathes’), see e.g. [56].
3. Tip frame
As in the non-resonant case, it is also possible to finda tip-based frame of reference:Φ = ϑ, X a = Z a Note that during the time interval T , the angle Φ needsto turn over ± π , whence ω = K Ω . Here, one can see that resonant meander is indeed aspecial case of Eqs. (8), in which the average value of ω tends to K Ω. In the tip frame, u ∗ represents a spiralthat performs no net rotation when Ψ is raised from 0 to2 π , it merely deforms around the same orientation. F. General case
The foregoing patterns can all be captured by notingthat each broken continuous symmetry delivers a collec-tive coordinate X M : x → X , y → X , θ → Φ , t → Ψ . We shall denote this as x M → X M ; there is one pairfor every broken continuous symmetry. We also convenethat { X M } = { X A , Φ , Ψ } , { X m } = { X a , Φ , Ψ } . The use of indices M vs. m thus distinguishes betweenthe moving and lab frames of reference, as did A vs. a for the translational modes only. Note that X a is tip po-sition in the lab frame, while x a is a lab frame Cartesiancoordinate that can label any point.The parameters of the solution x M are chosen from { x, y, θ, t } and X M are the collective coordinates. Thenwe find that all aforementioned sets of equations (5),(6), (7), and (8) that capture the evolution of particu-lar reaction-diffusion patterns can be written as U ( x a , t ) = u ∗ ( ρ A , Ψ ) ,∂ t X M = c M ( Ψ ) (10)given a time-dependent spatial coordinate transforma-tion ρ A ( x a , t ) = g ( x a , X M ( t )) that has as many degrees of freedom as broken Galileansymmetries by the solution.Note that the representation (10) is not unique. E.g.for a rigidly rotating spiral wave, one could choose U ( x a , t ) = u ∗ ( ρ A , Ψ ) that is periodic in its third argu-ment Ψ ≡ ϑ and ρ A identical to the lab frame coordi-nates x a . However, one can also choose to let u ∗ in theright-hand side of Eq. (10) only depend on ρ A and Ψ.For the actual calculations performed below, we workwith Eqs. (8) for a meandering spiral in a tip-based frameof reference, since they contain all other patterns and ref-erence frames as a special case, i.e. Eqs. (5), (6), (7)and(8). Those can be recovered from Eq. (8) by considering u ∗ independent of a collective coordinate. Thereafter,the evolution equation for that coordinate becomes irrel-evant and can be dropped. E.g. the theory of rigidlyrotating spirals follows from the meander case by taking u ∗ independent of Ψ and leaving out the evolution equa-tion for Ψ. From this case, the 1D traveling wave casefollows by considering u ∗ to be also independent of Φ and ρ , such that Eq. (5) follows.1 V. PROPERTIES OF THE LINEARIZEDPROBLEMA. Right-hand zero modes
Since the original RD Eq. (1) is invariant under spa-tial, rotational and temporal shifts, one has that, when asolution is rotated, or shifted in space or time, it is stilla solution to Eq. (1). This can be made explicit by sub-stituting u ( x + ˜ ε, y, t ) in Eq. (1), and similar for othercoordinates, yielding in first order in ˜ ε thatˆ (cid:96) ∂ U ∂x M = . (11)where ˆ (cid:96) = P ( ∂ x + ∂ y ) + F ′ ( U ) − ∂ t . (12)As expected, shifted values of the solution are right-handzero modes of the linearized operator ˆ (cid:96) associated to theRDE, in the lab frame. In several previous works, a moving frame was chosenin which the solution was either stationary [18, 32, 34,61, 65] or periodic [16]. This comes down to expressingthe derivative in Eq. (11) in the moving frame using thechain rule. From Eqs. (8) then follows ∂ t U = ∂ A u ∗ ∂ t ρ A + ∂ Ψ u ∗ ∂ t Ψ = ∂ A u ∗ [ ωR Aa ( Φ ) (cid:15) ab ( x b − X b ( t )) − R Aa ( Φ ) c a ]= Ω ∂ Ψ u ∗ − ω∂ θ u ∗ − c A ∂ A u ∗ . (13)Here, θ is the polar angle around the origin of the cho-sen moving reference frame. Thus, u ∗ obeys P ∆ u ∗ + F ( u ∗ ) + ω∂ θ u ∗ + c A ∂ A u ∗ − Ω ∂ Ψ u ∗ = . (14)Here, ω and c A may depend on Ψ (which is the case inthe tip frame of a meandering spiral wave). The linearoperator associated to Eq. (14) isˆ L = P ∆ + F ′ ( u ∗ ) + ω∂ θ + c A ∂ A − Ω ∂ Ψ Since the original RDE (1) is invariant in space, onehas that in an infinite medium, if U ( x a , t ) is a solution,so should U ( x a + ˜ ε a , t ) for any constant displacement ˜ ε a .We find in one spatial dimension thatˆ L ∂ ρ u ∗ = . In two spatial dimensions, we remark that ∂ a U = ∂ A u ∗ R Aa ( Φ ) , such that we find in first order in ˜ ε a :ˆ L ∂ A u ∗ = ω(cid:15) AB ∂ B u ∗ . (15)By taking complex combinations V ± = −( / )( ∂ x u ∗ ± i∂ y u ∗ ) , one obtains true eigenmodes of ˆ L , as ˆ L V ± =± iω V ± [14]. Note that in some of the chosen reference frames, ω may depend on Ψ, and the ‘eigenmodes’ of thesystem are 2 π -periodic functions of Ψ.For 2D patterns, we can state that a solution that isrotated around an angle ˜ ε around the current origin (i.e.spiral tip position or meander centre position) will stillbe a solution. With ∂ θ = (cid:15) AB ρ A ∂ B , one finds that, if U ( x a , t ) + ˜ ε∂ θ U ( x a , t ) is a solution, thenˆ L ∂ θ u ∗ = . We will denote − ∂ θ u ∗ as V Φ and − ∂ A u ∗ as V A .Thirdly, expressing that a time-shifted solution U ( x a , t + ˜ ε ) also solves (1) yieldsˆ L ∂ t u ∗ = . where ∂ t u ∗ , also denoted V Ψ , is given by Eq. (13).In the case of meandering spirals, V Φ is linearly in-dependent of V Ψ , since otherwise a shift in time wouldbe equivalent to simple rotation, and we would find our-selves in the non-meandering case.To summarise, equation (12) defines the set of N BS zero modes for the linearized operator ˆ (cid:96) in the laboratoryframe; in the comoving frame, according to (15), thisproduces a set of N BS eigenvaluesˆ LV ± = ± iω V ± , ˆ LV Φ = , ˆ LV Ψ = . In quantum field theories, bosons appearing due tospontaneous breakdown of continuous symmetries arecalled Goldstone bosons; extending the analogy to theclassical nonlinear field, the eigenfunctions correspond-ing to breakdown of continuous modes are sometimes re-ferred to as ‘Goldstone modes’.
B. Adjoint problem
Let us associate to ˆ L the operatorˆ L † = P T ∆ + F ′ T ( u ∗ ) − ω∂ θ − c A ∂ A + Ω ∂ Ψ . Note that u ∗ is 2 π -periodic in Ψ, and in the space offunctions 2 π -periodic in Ψ, ˆ L † is the adjoint operator toˆ L , in the sense that ⟪ ˆ L f ∣ g ⟫ = ⟪ f ∣ ˆ L † g ⟫ where ⟪ f ∣ g ⟫ = π ∫ π ⟨ f ∣ g ⟩ dΨ , (16a) ⟨ f ∣ g ⟩ = ∫ R N f H g d N x, (16b)where N = N = u ∗ is independent of Ψ, so the inner products in Eqs. (16a)and (16b) coincide, for any f and g defined by u ∗ .2Given that ˆ L † is the adjoint to ˆ L , we assume thatit also has N BS critical eigenmodes W M that are 2 π -periodic in Ψ:ˆ L † W ± = ∓ iω W ± , ˆ L † W Φ = , ˆ L † W Ψ = . (17)These ‘adjoint critical modes’ are known as sensitiv-ity functions or response functions (RFs) [19], as will beexplained in Sec. VI. C. Instant orthogonality of left and right criticalmodes
The set of critical adjoint modes can be normalised as ⟪ W M ∣ V N ⟫ = δ M N . Moreover, the orthogonality of critical and adjoint crit-ical eigenmodes holds instantaneously [16, 59, 66]: ⟨ W M ∣ V N ⟩ = δ M N , (18)Here, we show where the proof in [16] needs to be adaptedto accommodate for non-constant rotation rates ω ( Ψ ) and the case of resonance ( ∣ ω ∣ → Ω).Let us suppose that ˆ L † W M = λ M W M and ˆ L V N = λ N V N . Let us define the operatorsˆ L = ˆ L + Ω ∂ Ψ , ˆ L † = ˆ L † − Ω ∂ Ψ . (19)We note that ˆ L † is adjoint to ˆ L with respect to the innerproduct ⟨⋅ ∣ ⋅⟩ , which follows from integration by partsand the tempered nature of the adjoint eigenmodes W M .Denoting ⟨ W M ∣ V N ⟩ as I M N , it follows thatΩ ∂ Ψ I M N =⟨( ˆ L † − ˆ L † ) W M ∣ V N ⟩ + ⟨ W M ∣ ( ˆ L − ˆ L ) V N ⟩=( λ N − λ M ) I M N − ⟨ ˆ L † W M ∣ V N ⟩ + ⟨ W M ∣ ˆ LV N ⟩=( λ N − λ M ) I M N . (20)Hence I M N ( Ψ ) = I M N ( ) exp ⎛⎝ ∫ Ψ0 ( λ N ( z ) − λ M ( z )) d z ) Ω ⎞⎠ . (21)Now, if we have λ M = λ N , then I M N is constant. Then,setting ⟪ W M ∣ V N ⟫ = δ M N already yields (18). Else, if R ( λ M ) ≠ R ( λ N ) , and I M N ( ) ≠
0, then ∣ I M N ( Ψ )∣ growsor decays exponentially in time, which cannot happensince it should be 2 π -periodic. Therefore, I M N ( Ψ ) = R ( λ M ) = R ( λ N ) but I ( λ M ) ≠ I ( λ N ) , e.g. when considering the inner product between critical eigenmodes. In that case, λ M − λ N = iω ( Ψ ) (cid:96) with (cid:96) ∈ { , ± , ± } . Then, ∫ π ( λ N ( z ) − λ M ( z )) d z ) Ω = i(cid:96)β. (22)If not in the resonant case, β mod π ≠
0, and the ex-ponential factor in Eq. (21) cannot be periodic, whence I M N =
0. In the case of resonance, we are free to choosea non-rotating frame, where ω =
0, such that (18) stillholds. ◻ In the laboratory frame of reference, the critical modesare true zero modes, whence all λ m =
0, and the preser-vation of the inner product immediately follows from Eq.(20).
VI. SPATIOTEMPORAL DRIFT OF PATTERNSUNDER A SMALL PERTURBATIONA. Derivation of the drift equations
Using the ingredients defined above, it is possible topredict how a stable reaction-diffusion pattern (e.g. aplane wave or a spiral wave) reacts to a small perturba-tion. We mainly follow the derivation for rigidly rotatingspiral waves [36] but extend it to the case of meander. Incomparison to [16], we offer more flexibility in the frameof reference, such that also meandering spirals close toresonance can be treated.We start from the perturbed RDE: ∂ t u = P ( ∂ x + ∂ y ) u + F ( u ) + h ( x, y, t ) (23)where h = O( η ) is a small perturbation. A more genericform of perturbation can include dependence on the so-lution itself, i.e. h ( x, y, t, u , ∇) which however is reducedto (23) when the perturbation is evaluated at the unper-turbed solution, i.e. h ( x, y, t, u ∗ ( x, y, t ) , ∇) .If the initial state is close to a stable solution (i.e. trav-elling or spiral wave), the net effect of h will be to causea spatiotemporal drift of that pattern, which can be in-ferred from the collective coordinates. As before, we willpresent the result from the most general case of a mean-dering spiral wave (Eqs. (8)), from which all other casescan be inferred by eliminating some of the collective co-ordinates.Thus, we approximate the solution as u ( x a , t ) = u ∗ ( ρ A , Ψ ) + ˜ u ( x a , t ) , (24)where ˜ u = O( η ) . The coordinate transformation is nowgiven by ρ A = R Aa ( Φ ( t ))( x a − X a ( t )) , Ω τ = Ψ ( t ) . in which the collective coordinates’ temporal evolutionis perturbed by yet unknown drift terms v M , which are3also O( η ) : ∂ t X a = R Aa ( Φ )( c A + v A ) ,∂ t Φ = ω ( Ψ ) + v Φ ,∂ t Ψ = Ω + v Ψ . We can make the decomposition (24) unique by impos-ing that ⟨ W M ∣ ˜ u ⟩ = . (25)This is possible by shifting the solution. E.g. if oneapproximates at a given instance of time the solution u ∗ ( ρ , ρ , Ψ ) as u ∗ ( ρ + ˜ ε, ρ , Ψ ) + ˜ ε∂ u ∗ ( ρ , ρ , Ψ ) , then ⟨ W ∣ ˜ u ⟩ = ˜ ε ≠
0, and the proposed solution can be bettershifted to match the true solution, i.e. until Eq. (25)holds. A second interpretation of Eq. (25) is that thedeviation vector u should be orthogonal to the inertialmanifold [36].Next, plugging the ansatz (24) into the perturbed RDE(23) delivers: ∂ τ ˜ u − ˆ L ˜ u + v M V M = h + O( η ) . Note that by including the minus sign in the definitionof the GMs, i.e. V Φ = − ∂ θ u ∗ , V A = − ∂ A u ∗ , V Ψ = + ∂ Ψ u ∗ ,a plus sign appears in front of v M V M . Projection ontothe response function W M delivers, due to the meanderlemma (18): ⟨ W M ∣ ∂ τ − ˆ L ˜ u ⟩ − v M = ⟨ W M ∣ h ⟩ + O( η ) . The first term vanishes, due to the gauge condition (25)and the fact that ˆ L is adjoint to ˆ L with respect to ⟨⋅ ∣ ⋅⟩ : ⟨ W M ∣ ( ˆ L − ∂ τ ) ˜ u ⟩ = ⟨( ˆ L † + ∂ τ ) W M ∣ ˜ u ⟩= ⟨( ˆ L † + Ω ∂ Ψ ) W M ∣ ˜ u ⟩= ⟨ ˆ L W M ∣ ˜ u ⟩= λ M ⟨ W M ∣ ˜ u ⟩ = . Hence, we find the simple result: v M = ⟨ W M ∣ h ⟩ + O( η ) , (26)It turns out that at all times, the drift induced by a smallperturbation can be found by taking the overlap integralbetween the perturbation and the response function W M corresponding to that degree of freedom.This result justifies to call the critical adjoint eigen-modes ‘response functions’ [19]. In engineering terms,they are also the spatiotemporal ‘impulse response’ toa localised stimulus. This property can also be used toestimate RFs numerically or in future experiments [58].In the above-discussed moving frames of reference, theequations of motion become ∂ τ X A = c A ( Ψ ) + ⟨ W A ( Ψ ) ∣ h ⟩ ,∂ τ Φ = ω ( Ψ ) + ⟨ W Φ ( Ψ ) ∣ h ⟩ ,∂ τ Ψ = Ω + ⟨ W Ψ ( Ψ ) ∣ h ⟩ . which can also be summarised as [16]: ∂ t X M = c M ( Ψ ) + ⟨ W M ( Ψ ) ∣ h ⟩ Note the presence of a zeroth order motion c A ( Ψ ) ,which accounts for spiral tip motion even in the absenceof perturbations. Moreover, the rotation rate c ϑ ≡ ω isallowed to depend on the meander phase Ψ.To find the net drift effects, motion generally needsto be time-averaged in the laboratory frame of reference,starting from: ∂ t X a = R aA ( Φ ) c A ( Ψ ) + R aA ( Φ )⟨ W A ∣ h ⟩ ,∂ t Φ = ω ( Ψ ) + ⟨ W Φ ∣ h ⟩ ,∂ t Ψ = Ω + ⟨ W Ψ ∣ h ⟩ . (27)We will provide some elementary examples in Sec.VIII.At this stage, we can understand why different framesof reference may be used depending on the context. Letus suppose that we describe a process in which the per-turbation varies only slightly across the essential supportof the RFs of the system, say h (⃗ r ) = E a (⃗ r ) Q ∂ a u , where Q ∈ R n × n is a constant matrix acting in the space of re-active components (state variables). Then, the overlapintegrals will typically be approximated as [17, 18, 32, 35] ⟨ W M ∣ h ⟩ ≈ E a ( ⃗ X )⟨ W M ∣ Q ∂ a u ⟩ (28)and for non-homogeneous E N the next-to-leading orderwill be small only if the response functions W M are welllocalised near the point ⃗ X . Although the extent of theRF is fixed by the parameters of the model, the observercan thus describe the reduced system more accurately byan appropriate frame and tip choice. E.g. the interactionof a three-dimensional scroll wave in a detailed cardiacgeometry, the results using the tip frame are expectedto be more accurate than the center frame description.A similar argument holds for wave fronts, whose criticaladjoint eigenmodes are also localised around the wavefront. VII. GEOMETRIC INTERPRETATION:DYNAMICS ON THE INERTIAL MANIFOLD
The methods and results presented above can be rep-resented geometrically in the language of dynamical sys-tems on manifolds [55, 67]. The formal application of theidealized scheme presented below, as used in this paperas well as, explicitly or implicitly, in some other previousstudies, has well known technical difficulties related tothe facts that our phase space is a functional space, i.e.is infinite-dimensional, and that SE ( ) is a non-compactgroup. Some of the implications of these are discussedin the concluding section; for now, we proceed ignoringthese technical difficulties, in order to describe interest-ing physical phenomena, leaving rigorous mathematical4treatment for subsequent studies. The first part of thissection is an elaboration on the short discussion given inSec. II above. A. Phase space orbits
Every possible state of the system (1) at a fixed time,e.g. U (⃗ r ) , can be represented as a point of an infinite-dimensional phase space V . The set U (⃗ r, t ) obtainedby evolution from the initial condition U (⃗ r, t ) = U (⃗ r ) according to (1) is the orbit of U . For the formal devel-opment of the theory, we take t ∈ R , but we will not usebackward time evolution, which will be ill-posed for thesystems considered. Of the examples shown in Fig. 1,only rigidly rotating spiral waves and meandering spiralwaves with rational 2 π / β have closed orbits. B. Spatial vs. spatiotemporal symmetries
Since time is treated as an parameter in phase spacedynamics (i.e. it is not a ‘direction’ of the phase space), itwill be useful to distinguish between different symmetrygroups that include or not include time symmetry. As inSec. II, we take G = SE ( N )× R and H the largest contin-uous subgroup of G which leaves a given orbit invariant.The largest continuous subgroup of SE ( N ) leaving thegiven orbit invariant is denoted J . Then, we haveSE ( N ) ≤ G ≤ ≤ J ≤ H (29)The number of broken symmetries N BS , broken spacesymmetries N BSS and broken time symmetries N BT S bythe solution are then given by: N BS = dim ( G ) − dim ( H ) ,N BSS = dim ( SE ( N )) − dim ( J ) , (30) N BT S = N BS − N BSS . C. Quotient space
The group orbit of a state U (⃗ r ) under a group B isgiven by B U (⃗ r ) = { γ U (⃗ r ) ∣ γ ∈ B } . (31)It will be useful to study the dynamics of the system,disregarding spatial Euclidean symmetries. Therefore,we consider equivalence classes: U ∼ U iff there ex-ists γ ∈ SE ( N ) such that U = γ U . The space ofequivalence classes is called the quotient space, denoted Q = V/ SE ( N ) and contains as members the group orbitsunder SE ( N ) . The original dynamical system (1) induces a (reduced)dynamics on the states in Q .The equilibrium points of the dynamics in Q aretermed relative equilibria in V , as they correspond toequivalence classes under the action of the Galilean groupin V . Hence, from the examples discussed above, therigidly translating waves and rigidly rotating spirals arerelative equilibria of the reaction-diffusion system. (Forcomparison, the resting state is a ‘true’ equilibrium).The limit cycles of the dynamics in Q are termed rela-tive periodic orbits in V . Of the examples above, modu-lated traveling waves in one spatial dimension, biperiod-ically meandering spiral (both flower-like and star-like)and spirals in resonant meander are relative periodic or-bits.We will below further describe solution to the RD sys-tem that are relative equilibria or relative periodic orbits.These have an attractor A ⊂ Q of dimension 0 or 1. Theattractor is not necessarily unique to the medium. Forinstance, for the case of a 2D medium supporting spiralwaves, one already finds 4 different attractors: the rest-ing state, a traveling plane wave and spiral waves of twoopposite chiralities.
D. Inertial manifold
The presence of stable persisting patterns in the systemputs a special structure on the phase space. Considerthe set of states in V that belong to the same relativeequilibrium or relative periodic orbit: M = G A = { γ U (⃗ r ) ∣ U (⃗ r ) ∈ A , γ ∈ G } . (32)This set M (implicitly depending on a given choice for A )forms a manifold in the phase space V . From our previousassumption that ˆ (cid:96) has no eigenvalues with positive realpart, it follows that nearby trajectories are attracted withexponential speed to M , whence it is an inertial manifold [68] [69]. Since unidirectional attraction is implied here,the emergent structure of an inertial manifold is typicalfor dissipative systems.A given system (1) may have various inertial manifoldscorresponding to different prototypical solutions U , suchas the resting state, a propagating pulse, or a spiral wavewith either chirality.Inertial manifolds constitute a special case of slowmanifolds [67]. When the Euclidean symmetry is notbroken ( h = ), the dynamics on the inertial manifoldis trivial and given by the orbit of that solution. In thepresence of a perturbation, however, the solution will notfollow the original orbits anymore, but slowly drift fromone original orbit to the other. This regime has been thestarting point for many analytical studies of non-linearwave dynamics, e.g. [35, 65, 70–73], including extensionto 3D scroll waves [16–18, 32, 33, 74].A remarkable property is that in the high (infinite)-dimensional phase space, the inertial manifold associatedto the patterns above is finite-dimensional. Given that5the orbit U (⃗ r, t ) has dimension 1 and dim ( SE ( N )) = N ( N + )/
2, the maximal dimension of M would be N ( N + )/ +
1. However, a state U (⃗ r, t ) may be in-variant under a subgroup J of SE ( N ) , in which case J iscalled an isotropy subgroup . E.g. a plane wave in 2D isinvariant under translations along its wavefront.The elements of J correspond to the preserved sym-metries of the solution, such that the number of brokenspatial symmetries amounts to: N BSS = dim ( SE ( N )) − dim ( J ) . (33)For solutions that are relative equilibria, time evolu-tion can be represented by a Euclidean group action: U (⃗ r, t ) = γ U (⃗ r, ) with γ = γ ( t ) ∈ G for all t ∈ R . Inthis case, the time symmetry is not explicitly broken,and we say that N BT S =
0. For relative periodic orbits, U (⃗ r, t ) = γ ( t ) U (⃗ r, ) with γ ( t ) ∈ G only when t is aninteger multiple of the period T of the limit cycle in Q .In this case, we say that the number of broken time sym-metries N BT S = (M) = N BS = N BSS + N BT S . (34) E. Collective coordinates parameterize the inertialmanifold
By definition, the inertial manifold is flow-invariant, sothe dynamics on it can be described by a system of N BS coupled ordinary differential equations. Since it is a lo-cal attractor, this system will approximate the dynamicswithin the vicinity of that manifold.The different frames shown in this work serve the samepurpose: to uniquely characterize a dynamical state ofthe system. Different frame choices thus correspond todifferent parameterizations of the N BS -dimensional iner-tial manifold.Since every spatial collective coordinate is related to abroken symmetry, one has that ∂ U ∂X M = − ∂ U ∂ρ M , (35)where ρ M is a co-moving frame coordinate and the pro-totypical solution U ( x, t ) is assumed to depend on thecollective coordinates X M as parameters. E.g. since atraveling wave has the shape U ( x, t ) = u ∗ ( x − X ( t )) , onehas ∂ U ∂X = − ∂ u ∗ ∂x . (36)Similarly, for spiral waves ∂ U ∂X A = − ∂ u ∗ ∂ρ A , ∂ U ∂ϑ = − ∂ u ∗ ∂θ , (37)where θ is the polar angle in the chosen frame of referenceand ϑ is the spiral’s rotation phase. F. Tangent spaces
One can introduce tangent vectors to the inertial man-ifold by differentiating the state with respect to a collec-tive coordinate: V M = ∂∂X M U ( γ − ( X )⃗ r, t ) . (38)Here, γ − ( X ) has the effect of working in the quotientspace (or fully co-moving frame).If one differentiates with respect to the time-like col-lective coordinate (Ψ), one obtains a tangent vector tothe orbit of the state under time evolution according tothe RDE (1).These tangent vectors to the inertial manifold thus cor-respond to the Goldstone modes (GM) of the theory, orright-hand critical modes of the system, also called V M in sections V-VI. The non-critical modes of the systemthen point towards the other directions in phase-space,off the inertial manifold.In every point of the inertial manifold, one can alsointroduce a set of dual vectors W M , that correspond tothe response functions defined and used above. In everypoint of the inertial manifold it is possible to choose ⟨ W M ∣ V N ⟩ = δ M N . G. Effect of a small perturbation
Now, suppose one starts with a pattern state that ex-actly equals the prototypical solution U to the RDE.Hence the initial dynamics of this state is on the iner-tial manifold. Then, applying a small stimulus h at time t = h is carefully tuned, the state after applying theperturbation does not lie on the inertial manifold any-more. However, we assume that the wave solution U isdynamically stable, such that the resulting state evolvestowards the inertial manifold and converges to an orbitin it.Then, the resulting orbit will generally not be the samestate as would have been reached without the perturba-tion. That is, one retrieves the same solution (spiral ortraveling wave), but with different collective coordinates X M . A practical way of finding this shift, is to measurewhich part of the perturbation was tangential to the iner-tial manifold, and that is precisely what Eqs. (26) mean.The part of the perturbation that is orthogonal to theinertial manifold will decay over time, and only repre-sent transient dynamics (as long as the amplitude of theperturbation is small enough). H. Coordinate change in the tangent manifold
Geometrically speaking, the different coordinate sys-tems which we presented in Sections III and IV are alter-native ways to define collective coordinates, and therefore6alternative bases vectors for the tangent manifold. Withthat respect, we can state that ○ V M = ∂ ○ X N ∂ ˇ X M ˇ V N , (39) ○ W M = ∂ ○ X M ∂ ˇ X N ˇ W N . (40)such that ○ V = J ˇ V , (41)ˇ W = J H ○ W , (42)where J is the Jacobian of the coordinate transformationand J H is its Hermitian transpose. Eq. (42) is an exam-ple of the adjoint representation, and will be exemplifiedby Eqs. (50), (55). VIII. APPLICATION: SPIRAL DRIFT IN ACONSTANT EXTERNAL FIELDA. Motivation
A simple perturbation field that can be considered is aconstant vector field that couples to the gradient in oneor more state variables: h = ⃗ E ⋅ ⃗∇ Mu . (43)In chemical systems, this term can model an applied elec-trical field ⃗ E that induces drift of different ions u j withrespective mobilities M j that are found on the diagonal ofthe matrix M . The resulting phenomenon in a chemicalcontext is known as ‘electroforetic drift’.In 3D systems supporting spiral waves, the motion in-duced by a small filament curvature κ also generates aterm of the form (43), with ⃗ E = κ ⃗ N and M = P , where ⃗ N is the unit normal vector to the filament curve. B. Average drift of a circular core spiral: centerversus tip frame
For circular core spirals, the law of motion (27) withthe stimulus (43) yields ∂ t X a = R aA ( Φ ) c A + R aA ( Φ ) M AB R Bb ( Φ ) E b ∂ t Φ = ω + M Φ A R Aa ( Φ ) E a (44)with M M B = ⟨ W M ∣ MV B ⟩ . (45)The relations (44) hold both in the tip frame (with c A constant) and the center frame (with c A = t ( Φ ) instead of Φ ( t ) ,we immediately get t ( π ) − t ( ) = T + O( E ) (46)since the average of a rotation matrix over its angle van-ishes. Next, dividing the first equation of (44) by thesecond, we obtain ω∂ Φ X a = R aA ( Φ ) c A + (47) R aA ( Φ ) ( M AB − c A M Φ B ω ) R Bb ( Φ ) E b + O( E ) . Integrating over Φ on [ , π ] delivers the period-averageddrift velocity: V a = X a ( π ) − X a ( ) t ( π ) − t ( )= ( M AB − c A M Φ B ω ) ( δ ab δ BA + (cid:15) ab (cid:15) BA ) E b = γ E a + γ (cid:15) ab E b (48)with γ = ( M AA − c A M Φ A ω ) ,γ = (cid:15) BA ( M AB − c A M Φ B ω ) . (49)In the center frame, with c A =
0, this is the classical resultfrom [18]. However, does the calculation in the tip frameyield the same result?Here, we can make the discussion of Par. VII H ex-plicit. The center and tip frames are for a circular-corespiral related by ˇ ρ A = ○ ρ A − ζ A . Call the polar angle ineither frame ○ θ and ˇ θ ; the unperturbed spiral solution is ○ u ∗ ( ○ ρ A ) = ˇ u ∗ ( ˇ ρ A ) . Then, since ∂ ˇ θ = − (cid:15) BA ˇ ρ A ∂ B ˇ u ∗ , it fol-lows that ( ˇ V A ( ˇ ρ C ) ˇ V Φ ( ˇ ρ C ) ) = ( (cid:15) AB ζ B ) ⎛⎜⎜⎝ ○ V A ( ○ ρ C ) ○ V Φ ( ○ ρ C ) ⎞⎟⎟⎠ . (50)To find the RF transformation law, we consider a lo-calized perturbation of the j -th state variable at a givenposition: h = ( h k ) , h k = ηδ ( ρ A − ρ A ) δ ( t ) δ kj . From (27),we find in the center frame ∂ t ○ X a = ηR aA ( ○ ϑ ) δ ( t ) ○ W Aj ( ○ ρ A ) , (51) ∂ t ○ ϑ = ω + η ○ W Φ j ( ○ ρ A ) δ ( t ) , and in the tip frame ∂ t ˇ X a = R aA ( ˇ ϑ ) c A + ηR aA ( ˇ ϑ ) ˇ W Aj ( ˇ ρ A ) δ ( t ) , (52) ∂ t ˇΦ = ω + η ˇ W ϑj ( ˇ ρ A ) δ ( t ) . X a = ○ X a + R aA ( ○ ϑ ) ζ A with respect totime delivers ∂ t ˇ X a = ∂ t ○ X a − R aA ( ○ ϑ ) (cid:15) AB ∂ t ○ ϑ ζ B , (53) ∂ t ˇ ϑ = ∂ t ○ ϑ. After substituting Eqs. (51) and (52) into Eq. (53), weconclude from the zeroth order term in η that − (cid:15) AB ζ B ω = c A (54)and from the first order term in η that ⎛⎝ ˇ W A ( ˇ ρ C ) ˇ W Φ ( ˇ ρ C ) ⎞⎠ = ( − (cid:15) AB ζ B ) ⎛⎜⎝ ○ W A ( ○ ρ C ) ○ W Φ ( ○ ρ C ) ⎞⎟⎠ . (55)Hence, using Eq. (54) we can verify thatˇ γ = ⟨ ˇ W A − c A ω ˇ W Φ ∣ M ˇ V A ⟩ (56) = ⟨ ○ W A − (cid:15) AC ζ C ○ W Φ − c A ω ○ W Φ ∣ M ○ V A ⟩ = ○ γ and similar for γ . This result resolves an apparent para-dox: the expression for the filament tension γ is differentin the center frame and the tip frame, as shown by thefirst and second line of Eq. (56). However, both expres-sions are equal since the RFs and GMs appearing in theseexpression differ for both frames. The final result is thusindependent of the chosen frame of reference. C. Phase-locking of a meandering spiral under aconstant external field.
In systems where spiral waves are found where the an-gle β between subsequent petals of the meander flower isclose to a simple fraction of 2 π , the rotation phase maylock to the external field if that is strong enough [16].To capture this phenomenon quantitatively, we takethe minimally rotating finite-core frame of Fig. 3b.Then, when E = ∣ ⃗ E ∣ =
0, the X b ( t ) describes the mo-tion of a point moving on a circle with radius equal tothe mean radius of the meander flower. To recover theprecise tip motion from X b ( t ) one can take Z b ( t ) = X b ( t ) + R bB ( Φ ( t )) ζ B ( Ψ ( t )) . (57)Plugging Eq. (43) into Eq. (27) then delivers: ∂ t X a = R aA ( Φ ) c A + R aA ( Φ ) M AB ( Ψ ) R Bb ( Φ ) E b ,∂ t Φ = ω + M Φ B ( Ψ ) R Bb ( Φ ) E b ,∂ t Ψ = Ω + M Ψ B ( Ψ ) R Bb ( Φ ) E b . (58)The matrix element M M B is formally defined as inEq. (45), but now M ∈ { x, y, Φ , Ψ } . Note that functions M M B ( Ψ ) and R aA ( Φ ) are all 2 π -periodic. Due to ourchoice of coordinate system, c A and ω are non-zero butconstant.In our analysis, we first treat the case where β ≪ π ,i.e. near the 1 ∶ ω and E = ∣ ⃗ E ∣ are O( η ) , such that a locked solutionmay exist. Under those conditions, Ψ ( t ) will remain amonotonous function, and we can trade the time coordi-nate t for Ψ. Expanding up to linear order in η delivers:Ω ∂ Ψ X a = R aA ( Φ ) c A + O( η ) , (59a)Ω ∂ Ψ Φ = ω + M Φ B ( Ψ ) R Bb ( Φ ) E b + O( η ) . (59b)If there is phase-locking of the rotation angle Φ to themeandering phase Ψ, then equation (59b) will give a so-lution for Φ that remains in a bounded vicinity of somemean value Φ ∗ ,Φ ( Ψ ) = Φ ∗ + ˜Φ ( Ψ ) , ∣ ˜Φ ( Ψ )∣ < const . We will be looking for the simplest case of 1:1 phaselocking, such that Φ ( Ψ + π ) ≡ Φ ( Ψ ) . From Eq. (59b), the rate-of-change of Φ is also small(i.e. O( η ) ), while the rate-of-change of Ψ is Ω = O( ) .Therefore, we can assume that ˜Φ = O( η ) during onemeander cycle, and average over the cycle. DenotingΦ ( πn ) as Φ n , we obtain the difference equation:Φ n + − Φ n T = ω + M Φ B R Bb ( Φ n ) E b . (60)Then, a phase-locked angle is a fixpoint of Eq. (60).Without loss of generality, we can take ⃗ E = E ⃗ e x here, toshow that this condition can be formulated as: ω + AE cos ( Φ ∗ − α ) = Ae − iα = M Φ1 + i M Φ2 . A necessary condition for a phase-locked solution toexist is then ∣ EA / ω ∣ >
1, i.e. [16] ∣ E ∣ > E crit = ωA + O( η ) . From Eq. (61), the rotation angle around which thelocked equation will equilibrate isΦ ∗ = α + arccos ( ωAE ) . The second solution to Eq. (61), namely α − arccos ( ω / AE ) , is unstable.Since in the phase-locked state, Φ = Φ ∗ + O( η ) , Eq.(59a) can be readily integrated to find the net spatialdisplacement during a phase-locked meander cycle: X a ( π ) − X a ( ) = R aA ( Φ ∗ ) c A T + O( η ) , (62)8 FIG. 5. Magnitude of spiral wave drift velocity in Barkley’smodel ( b = . ς = . a = .
58, black),and close to resonance ( a = .
63, red). Close to the phase-locking threshold, the drift speed approaches the orbital ve-locity c = ωr , with r the average radius of the meander flower. where T = π / Ω. From ∂ t Ψ = Ω + O( η ) , one finds more-over that the time needed to make a full meander cycleis T + O( η ) , such that the net drift speed during phaselocking is V a = X a ( π ) − X a ( ) T + O( η )= R aA ( Φ ∗ ) c A + O( η ) . (63)In words, it can be concluded that the drift speed of ameandering spiral phase-locked to a constant vector field ⃗ E is of the order of the ‘orbital’ velocity c , i.e. the netspeed at which the tip traverses the rim of the meanderflower. Importantly, within the locking region the driftvelocity is not proportional to the perturbation strength E ! In terms of the velocity parallel or perpendicular tothe applied field, Eq. (63) reads: V ∥ = c cos Φ ∗ + O( η ) , V ⊥ = − c sin Φ ∗ + O( η ) , a result which was already quoted in [16]. The magnitudeof drift velocity is therefore expected to be V = c + O( η ) .Fig. 5 confirms that this is the case for simulations inBarkley’s model.Until here, we have described phase-locking of a singlemeander cycle to an external field, which happens when ∣ ω ∣ ≪ Ω. Since Ω = π / T and ω = β / T with β the anglebetween subsequently visited petal tips, phase-locking isfound generally when β ≈ πp, p ∈ Z . In some cases, a higher-order resonance may occur, i.e. β ≈ πp / q, p, q ∈ Z , q ≠ p, q ).In this case, q meander cycles can be re-interpreted as asingle cycle to which the theory (27)-(63) can be applied,with ω replaced by qω − p Ω. D. Regime close to phase-locking
It could happen that the applied field ⃗ E is not strongenough to enforce phase-locking, which occurs when ∣ E ∣ < E crit . In that case, the net drift can be found from theinstant equations of motion (27) by averaging of one sortor another. In the case of circular core spirals, this wasdone by sliding time averaging [36] or by imposing condi-tion of periodicity on the perturbed solution that in turnrequiring a solvability condition [20]. Since for meander-ing spirals the unperturbed solution is not periodic butbiperiodic, the situation here requires a bit more care.In [58] we considered the case below but not too closeto phase-locking threshold, and used Fourier series. Forthe case ∣ E ∣ ≲ E crit , this does not work, so we resort tothe more robust method, by sliding averages. First wedo averaging over the meander phase:Φ ( Ψ ) = π ∫ Ψ + π Ψ − π Φ ( Ψ ′ ) dΨ ′ . Then, if ∣ ω ∣ ≪ Ω, ∂ t Φ = ω + M Φ A R Aa ( Φ ) E a + O( E ) . (64)As above, we can write the second term in the right-handside as EA cos ( Φ − α ) . Let us denote the time needed tocomplete a full meander flower as T ∗ = πω = πβ T . Under a perturbation this value changes to T . Then,straightforward integration of Eq. (64) delivers: T = T ∗ √ − ( E / E crit ) . Therefore, if E approaches E crit from below, the timeneeded to complete a full meander flower diverges. Be-yond E crit the spiral fails to complete a meander cycleand is phase-locked. E. Average drift speed of a non-phase-lockedmeandering spiral
We continue our analysis of Eq. (27) in the regimewhere phase-locking does not occur, i.e. when β / π is9 FIG. 6. Non phase-locked drift of a meandering spiral inBarkley’s model due to a constant external field ⃗ E = E ⃗ e x ,for different field strengths: (a) E = .
01, (b) E = .
02, (c) E = . Q = P . Paneld) shows the parallel and perpendicular components of drift(data points), which closely match the theoretical predictionsgiven by Eq. (71) (lines).FIG. 7. Same as Fig. 6, but for the Fenton-Karma guinea pigcardiac model from Fig. 1, which has a linear core. Appliedfield strengths are (a) E = .
005 mm − , (b) E = .
025 mm − ,(c) E = .
06 mm − . not a simple fraction or when E is sufficiently small, say E = O( η ) .For any of the reference frames shown in Fig. 3, wecan in a first step integrate Eqs. (58) to find the netdisplacement in space, time and rotation angle duringthe single meander cycle, say the n -th. At the start ofthis cycle, taking place at t = t n the collective coordinatesof the spiral are denoted X an , Φ n , Ψ n .For t ∈ [ t n , t n + ] , we can write:Ψ ( t ) = Ψ n + Ψ ∗ ( t ) + ˜Ψ ( t ) Φ ( t ) = Φ n + Φ ∗ ( t ) + ˜Φ ( t ) X a ( t ) = X an + X a ∗ ( t ) + ˜ X a ( t ) . In these right-hand sides, the first term is the initialvalue at the start of the n -th cycle, the second term isthe evolution in the unperturbed case, and the third termthe correction to it, of O( E ) . Without loss of generality,we can set the n -th cycle to start at Ψ = Ψ n = πn .The unperturbed evolution is given byΨ ∗ ( t ) = Ω ( t − t n ) Φ ∗ ( t ) = ∫ tt n ω ( Ψ ∗ ( t ′ )) d t ′ X a ∗ ( t ) = ∫ tt n R aA ( Φ ∗ ( t ′ )) c A ( Ψ ∗ ( t ′ )) d t ′ Due to the group structure of rotations, we can write R Aa ( Φ ) = R Aα ( Φ ∗ ) R αβ ( ˜Φ ) R βa ( Φ n ) ,R αβ ( ˜Φ ) = δ αβ cos ˜Φ + (cid:15) αβ sin ˜Φ = δ αβ + (cid:15) αβ ˜Φ + O( η ) . (65)We will below useM µα ( Ψ ∗ ( t )) = M µA ( Ψ ∗ ( t )) R Aα ( Φ ∗ ( t )) , µ ∈ { Φ , Ψ } M αβ ( Ψ ∗ ( t )) = R αA ( Φ ∗ ( t )) M AB ( Ψ ∗ ( t )) R Aα ( Φ ∗ ( t )) ,c α ( Ψ ∗ ( t )) = R αA ( Φ ∗ ( t )) c A ( Ψ ) (66)Note that these functions are not periodic in Ψ any-more.With Eqs. (65)-(66), the evolution equation for Ψ be-comes ∂ t Ψ = Ω + M Ψ β ( Ψ ∗ ( t )) R βb ( Φ n ) E b + O( η ) . (67)Since Φ n is constant, Eq. (67) yields the duration ofthe meander period in the presence of a constant externalfield ⃗ E : T n = t n + − t n = ∫ t n + t n d t = ∫ π dΨΩ + M Ψ α R αa ( Φ n ) E a = T ( − Ω − M Ψ α R αa ( Φ n ) E a ) + O( η ) . (68)Note that the resulting period depends on Φ n , i.e. onthe orientation of the trajectory relative to the appliedfield ⃗ E .0The rotation phase during one cycle can be found byintegrating ∂ Ψ Φ = ∂ t Φ / ∂ t Ψ:Φ ( Ψ ) − Φ n = ∫ ΨΨ n ω ( Ψ ) + M Φ α ( Ψ ) R αa ( Φ n ) E a Ω + M Ψ α ( Ψ ) R αa ( Φ n ) E a dΨ = Φ ∗ + N Φ α R αa ( Φ n ) E a + O( η ) , where N Φ α ( Ψ ) = Ω − ∫ ΨΨ n m α ( Ψ ′ ) dΨ ′ ,m α ( Ψ ) = M Φ α ( Ψ ) − Ω − T ω M Ψ α ( Ψ ) . Evaluation at Ψ n + delivers the change in rotation angleover one meander period: β n = Φ n + − Φ n = β + T m α R αa ( Φ n ) E a + O( η ) , (69)Hence, an applied field will change the angle betweenconsecutively visited petals of the meander flower.Finally, we can find the net spatial displacement duringa meander cycle: d an = X an + − X an = d a + T R aα ( Φ n ) m αβ R βb ( Φ n ) E b + O( η ) , wherem αβ = M αβ − Ω − [ c α M Ψ β − (cid:15) αγ c γ N Φ β ] . (70)Thus, the net drift of a meandering spiral in an exter-nal field depends on the angle Φ between the meanderpattern and the applied field ⃗ E .The average drift speed, measured over one meanderperiod (i.e. one petal of the flower) will be: V an = d an / T n = R aα ( Φ n ) Q αβ c α + R aα ( Φ n ) R βb ( Φ n ) E b + O( η ) . with Q αβ = m αβ + c α M Ψ β ΩIf we are outside the phase-locking regime, the ratio β / π = ω / Ω will not be close to a fraction p / q with p and q small integers. Therefore, after several rotations, the an-gle Φ n will have taken different values that are uniformlydistributed over the unit circle. In that approximation,the time average of a quantity (net rotation, displace-ment or time lapse) will be that quantity averaged overall possible phases Φ n at the start of the meander cycle: ⟨ T ⟩ = lim m →∞ t n + m − t n m = π ∫ π T n ( Φ n ) dΦ n , ⟨ β ⟩ = lim m →∞ Φ n + m − Φ n m = π ∫ π β n ( Φ n ) dΦ n , ⟨ d a ⟩ = lim m →∞ X an + m − X an m = π ∫ π d an ( Φ n ) dΦ n , ⟨ V a ⟩ = lim m →∞ X an + m − X an t n + m − t n = π ∫ π V an ( Φ n ) dΦ n . Inserting Eqs. (68), (69), the average duration androtation angle over a meander cycle satisfy ⟨ T ⟩ = T + O( η ) , ⟨ β ⟩ = β + O( η ) . For the average drift speed, we then find ⟨ V a ⟩ = Γ E a + Γ (cid:15) ab E b , (71)where Γ =
12 Q αα , (72)Γ = (cid:15) αβ Q βα , Since the trace of a matrix is invariant and (cid:15) = ( (cid:15) αβ ) is the generator of rotations, we furthermore find that forany 2 × P and rotation matrix R holdsTr ( RPR T ) = Tr ( R T RP ) = Tr ( P ) , Tr ( (cid:15) RPR T ) = Tr ( R T (cid:15) RP ) = Tr ( R T R (cid:15) P ) = Tr ( (cid:15) P ) . Hence Q αα = Q AA ,(cid:15) αβ Q βα = (cid:15) AB Q BA . Expression (71) was also computed in the center framein [16] using a double Fourier series. We find the sameresult here using the discrete maps for X an and Φ n . FromFigs. 6 and 7 it can be seen that the net drift velocity iswell described by Eq. (71) with coefficients (72). IX. DISCUSSION
In this work, we have introduced different coordinatesystems that allow to describe the dynamics of one- andtwo-dimensional wave patterns. The origin of the chosencoordinate system is interpreted as the wave front or spi-ral wave tip position, or as the filament position for 3Dscroll waves.For spiral waves, different coordinate systems producethe same result where they are applicable, but they candiffer in their applicability areas. So, the “minimally ro-tating” frame of reference is not good for describing driftof meandering spirals for parameter values close to theequivariant Hopf bifurction (Winfree’s ∂M boundary),and one needs to use the “co-rotating finite-core” frameinstead. Similarly, the “center” frame is inadequate forthe parameters in the vicinity of the “resonant meander”parameters, and one is much better with the “tip” framethere.From the simple application of electroforetic drift al-ready, it can be seen that the time-scale at which themotion is analysed is important. Here, we have demon-strated that different time-averaging strategies can be1taken in addition to the Fourier approach of [16]: analysisof the return map and a sliding average are possible.We have deliberately included a brief description ofthe geometrical viewpoint on phase space dynamics inSec. VII, which shows that our technical computationsare part of a simpler geometric theory.We have given two applications of the response func-tion framework here, both related to spiral wave driftin a constant external field. First, we have shown thatthe mean drift speed is equal in the center frame andfinite-core frame, as both overlap integrals produce thesame physical tension coefficients. Secondly, we havecomputed that for phase-locked meander, the drift speedis not proportional to the applied field strength, but closeto the orbital velocity along the meander flower.As always, the limitations of the presented work sug-gest the directions for future research. We mentionedin the beginning of Section VII that the mathematicalfoundation of our formal approach has known outstand-ing technical difficulties, e.g. to describe spiral waves thereaction-diffusion system has to be considered in a Ba-nach space in which the action of the Euclidean group isnot even continuous, and that the effective spatial local-ization of the response functions is only an empirical fact.A possible route through the first obstacle was shownby [54] which informally can be summarized thus: ratherthan considering the whole Banach space, one ought tofocus on its part occupied by actual spiral-wave solutions, and in that part, the group acts “nicely”. For the sec-ond difficulty, a promising rigorous result was obtainedin [75]: that response functions of one-dimensional ana-logues of spiral waves are indeed exponentially localizedin space.From the physical viewpoint, more interesting are thelimitations related to the types of perturbations consid-ered. Of special interest is wave propagation in confinedirregular geometries such as cardiac tissue. Previousworks have shown that these can also to some extent bedescribed by a perturbative approach [15, 76], althoughthe boundary perturbation may take the shape of a Diracdistribution and is therefore not small at all.Note that phase-locking of meandering spirals causedby perturbations that are purely time-dependent havebeen considered before [77, 78]. Our approach outlinedin [16] and detailed in here is potentially by far more uni-versal. The simple examples considered so far are admit-tedly only first steps, and consideration of more realisticproblems, such as cardiac tissue with its anisotropy andspatial heterogeneity, promises many new discoveries.This research was supported in part the EPSRCgrants EP/E018548/1, EP/D074789/1, EP/P008690/1,EP/N014391/1 and EP/E016391/1 (UK) and NationalScience Foundation Grants No. NSF PHY-1748958, NIHGrant No. R25GM067110, and the Gordon and BettyMoore Foundation Grant No. 2919.01 (USA). [1] W. Jahnke, C. Henze, and A. Winfree, Nature , 662(1988).[2] S. Nettesheim, A. von Oertzen, H. Rotermund, andG. Ertl, J. Chem. Phys. , 9977 (1993).[3] J. Lechleiter, S. Girard, E. Peraltal, and D. Clapham,Science , 123 (1991).[4] F. Siegert and C. Weijer, Proc. Natl. Acad. Sci. USA. ,6433 (1992).[5] N. Gorelova and J. Bures, J. Neurobiol. , 353 (1983).[6] M. Allessie, F. Bonke, and F. Schopman, Circ. Res. ,54 (1973).[7] R. Gray, A. Pertsov, and J. Jalife, Nature , 75 (1998).[8] F. Witkowsky, L. Leon, P. Penkoske, W. Giles, M. Spano,W. Ditto, and A. Winfree, Nature , 78 (1998).[9] M. Haissaguerre, M. Hocini, A. Denis, A. J. Shah, Y. Ko-matsu, S. Yamashita, M. Daly, S. Amraoui, S. Zeller-hoff, M.-Q. Picat, A. Quotb, L. Jesel, H. Lim, S. Ploux,P. Bordachar, G. Attuel, V. Meillet, P. Ritter, N. Derval,F. Sacher, O. Bernus, H. Cochet, P. Jais, and R. Dubois,Circulation , 530 (2014).[10] S. M. Narayan, D. E. Krummen, K. Shivkumar, P. Clop-ton, W.-J. Rappel, and J. M. Miller, Journal of theAmerican College of Cardiology , 628 (2012).[11] I. Biktasheva, Y. Elkin, and V. Biktashev, Phys. Rev. E , 2656 (1998).[12] I. V. Biktasheva and V. N. Biktashev, J. Nonlin. Math.Phys. , 8 (2001).[13] V. Hakim and H. Henry, Phys. Rev. E , 046235 (2002).[14] I. Biktasheva, D. Barkley, V. Biktashev, G. Bordyugov, and A. Foulkes, Phys. Rev. E , 056702 (2009).[15] C. D. Marcotte and R. O. Grigoriev, Chaos , 063116(2015).[16] H. Dierckx, H. Verschelde, A. Panfilov, I. Biktasheva,and V. Biktashev, Phys. Rev. Lett. , 258101 1 (2017).[17] J. Keener, Physica D , 269 (1988).[18] V. Biktashev, A. Holden, and H. Zhang, Phil. Trans. R.Soc. Lond. A , 611 (1994).[19] I. Biktasheva and V. Biktashev, Phys. Rev. E , 026221(2003).[20] I. V. Biktasheva, D. Barkley, V. N. Biktashev, and A. J.Foulkes, Physical Review E , 066202 (2010).[21] V. Biktashev and A. Holden, Physica D , 342 (1998).[22] M. Nicol, I. Melbourne, and P. Ashwin, Nonlinearity ,275 (2001).[23] P. Ashwin, I. Melbourne, and M. Nicol, Physica D ,364 (2001).[24] R. Clayton, O. Bernus, E. Cherry, H. Dierckx, F. Fenton,L. Mirabella, A. Panfilov, F. Sachse, G. Seemann, andH. Zhang, Prog Biophys Molec Biol , 22 (2011).[25] V. Biktashev and M. Tsyganov, Phys. Rev. Lett. ,134101 1 (2011).[26] R. FitzHugh, Biophys. J. , 445 (1961).[27] D. Barkley, Physica D , 61 (1991).[28] F. Fenton and A. Karma, Chaos , 20 (1998).[29] V. Zykov, A. Mikhailov, and S. M¨uller, Phys. Rev. Lett. , 3398 (1997).[30] R. Clayton, E. Zhuchkova, and A. Panfilov, Prog Bio-phys Molec Biol , 378 (2005). [31] R. A. Gray, J. Wikswo, and N. Otani, Chaos , 033118(2009).[32] H. Verschelde, H. Dierckx, and O. Bernus, Phys. Rev.Lett. , 168104 (2007).[33] V. Biktashev, D. Barkley, and I. Biktasheva, Phys. Rev.Lett. , 058302 (2010).[34] H. Dierckx, O. Selsil, H. Verschelde, and V. Biktashev,Phys Rev Lett , 174102 (2012).[35] H. Dierckx, E. Brisard, H. Verschelde, and A. V. Pan-filov, Physical Review E , 012908 (2013).[36] V. Biktashev and A. Holden, Chaos, Solitons and Frac-tals , 575 (1995).[37] A. Panfilov and H. Dierckx, “Theory of cardiac arrhyth-mias,” in Cardiac electrophysiology. From cell to bedside,7nth edition , edited by D. Zipes, J. Jalife, and W. Steven-son (Elsevier, London, 2017) pp. 325–335.[38] A. Winfree, Science , 937 (1973).[39] O. Steinbock, V. Zykov, and S. M¨uller, Nature , 322(1993).[40] M. Yamazaki, S. Mironov, C. Taravant, J. Brec, L. M.Vaquero, K. Bandaru, U. Avula, H. Honjo, I. Kodama,O. Berenfeld, and J. Kalifa, Cardiovasc. Res. , 48(2012).[41] I. Efimov, V. I. Krinsky, and J. Jalife, Chaos, Solitonsand Fractals , 513 (1995).[42] M. Courtemanche, R. Ramirez, and S. Nattel,Am.J.Physiol. , H301 (1998).[43] V. Biktashev and A. Holden, Proc Roy Soc London B , 1373 (1996).[44] A. Bueno-Orovio, E. Cherry, and F. Fenton, J TheorBiol , 544 (2008).[45] O. E. R¨ossler and C. Kahlert, Z. Naturforsch. , 565(1979).[46] V. S. Zykov, Biofizika , 862 (1986).[47] A. Winfree, Chaos , 303 (1990).[48] D. Barkley, M. Kness, and L. Tuckerman, Phys Rev A , 2489 (1990).[49] A. Karma, Physical Review Letters , 2824 (1990).[50] D. Barkley, Phys.Rev.Lett. , 164 (1994).[51] C. Wulff, Theory of Meandering and Drifting SpiralWaves in Reaction-Diffusion Systems , Ph.D. thesis, FreieUniversit¨at Berlin (1996).[52] B. Fiedler, B. Sandstede, A. Scheel, and C. Wulff, Doc-umenta Mathematica , 479 (1996).[53] M. Golubitsky, V. LeBlanc, and I. Melbourne, Nonlin.Sci. (1997).[54] B. Sandstede, A. Scheel, and C. Wulff, J. DifferentialEquations , 122 (1997).[55] V. Biktashev, A. Holden, and E. Nikolaev, Int J BifurcChaos , 2433 (1996).[56] A. Foulkes and V. Biktashev, Phys Rev E , 046702(2010).[57] H. Dierckx and H. Verschelde, Phys. Rev. E , 062907(2013).[58] H. Dierckx, H. Verschelde, and A. Panfilov, Chaos ,093912 1 (2017).[59] C. D. Marcotte and R. O. Grigoriev, Chaos , 093107(2016).[60] N. Manz and O. Steinbock, Chaos , 037112 1 (2006).[61] H. Dierckx, O. Bernus, and H. Verschelde, Phys D ,941 (2009).[62] B.-W. Li, M.-C. Cai, H. Zhang, A. Panfilov, and H. Dier-ckx, J. Chem. Phys. , 184901 1 (2014).[63] V. Zykov, Simulation of wave processes in excitable media (Manchester University Press, Manchester, 1987).[64] C. Wulff,
Theory of Meandering and Drifting SpiralWaves in Reaction-Diffusion Systems , Ph.D. thesis, FreieUniversit¨at Berlin (1996).[65] J. Keener, SIAM J Appl Math , 1039 (1986).[66] A. Foulkes, Drift And Meander Of Spiral Waves , Ph.D.thesis, University of Liverpool (2009).[67] A. Debussche and R. Temam, Applied mathematics let-ters , 73 (1991).[68] R. Temam, Infinite-Dimensional Dynamical Systems inMechanics and Physics (Springer, New York, 1988).[69] Technically speaking, this is not precisely true: thereare continuous branches of the spectrum that reach theimaginary axis, as a result of which the attraction is notexponential. However, on the “physical level of rigour”,these continuous branches correspond to large-scale per-turbations of the periphery of the spiral wave, so if weconsider a finite vicinity of the spiral core, the approachis still sound. This is one of the outstanding technicaldifficulties we mentioned earlier.[70] K. Kuramoto, Prog. Theor. Phys. , 1885 (1980).[71] H. Henry, Phys Rev E , 026204 (2004).[72] H. Dierckx, O. Bernus, and H. Verschelde, Phys RevLett , 108101 (2011).[73] J. Loeber and H. Engel, Phys. Rev. Lett. , 148305(2014).[74] H. Henry and V. Hakim, Phys Rev E , 046235 (2002).[75] B. Sandstede and A. Scheel, SIAM J. Applied DynamicalSystems , 1 (2004).[76] I. Biktasheva, H. Dierckx, and V. Biktashev, Phys. Rev.Lett. , 068302 (2015).[77] S. Grill, V. Zykov, and S. M¨uller, J. Phys. Chem. ,19082 (1996).[78] R.-M. Mantel and D. Barkley, Physical Review E ,4791 (1996). APPENDIX: LIST OF MAIN NOTATIONS
Note that some symbols denote different objects in dif-ferent places; we hope their meaning is sufficiently clearfrom the context. f ∗ unperturbed value of f ˜ f perturbed value of f ○ f value of f in the center frameˇ f value of f in the tip frame f average of f over a cycle f sliding average of f ⟨ . ∣ . ⟩ inner product in space ⟪ . ∣ . ⟫ cycle-averaged inner product ⟨ . ⟩ average over uniformly distributed initial con-ditions a, b, ... spatial indices in the lab frame of reference,in { } A, B, ... spatial indices in the moving (rotating) frameof reference a, b, ς kinetic parameters of the Barkley model3 α, β spatial indices in a frame that follows the un-perturbed dynamics α, β, γ kinetic parameters of the FitzHugh-Nagumomodel β angle between petals of the meander flower c
1D frame velocity c A velocity components of the moving frame, ex-pressed in the moving frame C middle of unperturbed spiral wave tip trajec-tory∆ Laplacian operator d a net displacement by the spiral wave tip overone meander period (cid:15) Levi-Civita symbol˜ ε a components of infinitesimal translation vector ⃗ E external field η order of magnitude of the external perturba-tion F column-vector of n nonlinear functions repre-senting reaction kineticsΦ rotation angle of the moving frame γ element of transformation group γ scalar filament tension (circular core case) γ pseudoscalar filament tension (circular corecase) G system’s symmetry groupΓ scalar filament tension (meander case)Γ pseudo-scalar filament tension (meandercase) H solutions’s spatiotemporal symmetry group h perturbation to the reaction-diffusion system j, k spatial indices in the lab frame of reference J Jacobian matrix for the collective coordinates J solutions’s spatial symmetry group K spiral wave chirality (-1 or +1)ˆ (cid:96) linearized operator in the lab frameˆ L linearized operator in the moving frame m, n, ... generalized indices in the lab frame of refer-ence, in { x, y, Φ , Ψ } M , N , ... generalized indices in the moving frame, in { , , Φ , Ψ }M coupling matrix n number of state variables in the reaction-diffusion system N number of spatial dimensions N BS number of broken continuous spatiotemporalsymmetries N BSS number of broken continuous spatial symme-tries N BT S number of broken continuous temporal sym-metries (0 or 1) O origin of the lab frame ω instantaneous rotation velocity of the spiralwave ω cycle-averaged rotation velocity of the spiralwaveΩ frequency of temporal cycle completion (Ω = π / T ) P matrix of diffusivitiesΨ temporal phase of the relatively periodic pat-tern / meander phase Q quotient space ⃗ r position vector in the lab frame R aA rotation matrix ρ spatial coordinate in the moving frame of apoint ( ρ = x − X ) t time in the lab frame τ time in the moving frame T time needed to describe a full meander flower T duration of a meander cycle / relative orbit ϑ polar angle in the moving frame ( ϑ = θ − Φ ) θ polar angle in the lab frame u column-vector of state variables u j j -th state variable u ∗ j pinning constant u ∗ unperturbed spiral solution in the movingframe U solution to the reaction-diffusion equation inthe lab frame (column-vector) V M critical eigenmode V phase space v M instantaneous drift speed V a components of average spiral wave drift speed w m response function / critical adjoint eigenmodein lab frame W M response function / critical adjoint eigenmodein a moving frame x a lab frame coordinates of a point x, y lab frame coordinates of a point X, X a lab frame coordinate(s) of the origin of themoving frame X lab frame coordinate of the origin of the mov-ing frame X M collective coordinates X origin of the moving frame χ polar angle in kinetics phase space Z wave front / spiral wave tip position Z wave front / spiral wave tip coordinate in thelab frame ζζ