Measurement and structure of spiral wave response functions

The rotating spiral waves that emerge in diverse natural and man-made systems typically exhibit a particle-like behaviour since their adjoint critical eigenmodes (response functions) are often seen to be localised around the spiral core. We present a simple method to numerically compute response functions for circular-core and meandering spirals by recording their drift response to many elementary perturbations. Although our method is computationally more expensive than solving the adjoint system, our technique is fully parallellisable, does not suffer from memory limitations and can be applied to experiments. For a cardiac tissue model with the linear spiral core, we find that the response functions are localised near the turning points of the trajectory.


I. INTRODUCTION
Spiral waves are a prime example of self-organisation and have been observed in diverse systems, including catalytic oxidation 1 and oscillating chemical reactions. 2Moreover, they seem to be ubiquitous in biological systems, where they organise rhythmic activity from the macroscopic scale (e.g., amoeba morphogenesis 3 ) down to intracellular kinetics (e.g., intracellular calcium waves 4 ).Interestingly, at the tissue scale, they have been observed in neural 5 and cardiac tissue, [6][7][8] where they are thought to underlie cardiac arrhythmias.
The observation that the electrical activity in the heart during arrhythmias can be organised as spiral waves has motivated their study for decades, and in the medical community, they are known as "rotors."Recent clinical results suggest that ablating the tissue at rotor locations may cure fibrillation of the cardiac atria. 9However, if the rotors occur in the cardiac ventricles, the pumping of blood is quickly lost, and an electrical defibrillating shock needs to be administered within minutes.These examples show that a better knowledge of cardiac rotors and their dynamics may lead to better treatment or defibrillation techniques.
It was noted in some of the systems above and numerical simulations thereof that spiral waves can be remarkably stable structures.One explanation to this observation is that spiral waves are "topologically protected" structures, in the sense that when the activation phase is defined, they exhibit a phase singularity close to their rotation center, i.e., in their core region. 8,10Another possible explanation is that the tail of the spiral wave sweeps the surrounding space, resetting the medium properties to its resting values when it is passed.This way, external disturbances are strongly attenuated before they can affect the dynamics at the spiral wave core.
When the spiral wave is a dynamical attractor of the system in a two-dimensional planar geometry, its instantaneous state can be well approximated by a small set of collective coordinates, which tell how the Euclidean symmetries of the system are broken by the spiral wave solution.For a rigidly rotating spiral, we can track the spiral wave tip by computing the intersection of two isolines of variables. 11Thereby, one obtains the spiral wave tip position ðXðtÞ; YðtÞÞ in a Cartesian frame and its rotation phase UðtÞ.If the spiral wave solution is quasi-periodic (meandering), the solution is only periodic modulo a Euclidean transformation.For non-resonant meander, this Euclidean transformation can be taken to be a rotation over the rotation phase / around the centre of the meander flower.In the new frame of reference (i.e., in the quotient system of the dynamical system), the spiral solution becomes time-periodic, which can be converted to a meander phase WðtÞ ¼ X 0 t þ W 0 , such that the solution becomes 2p-periodic in W. The meander phase WðtÞ tells how far the solution has gone through the meander period.Below, we will group the collective coordinates as X l ¼ ðX; Y; U; WÞ and W is only included for meandering spirals.
For definiteness, we work with spiral wave solutions uðr; tÞ : X Â R !R N , i.e., the state u at every point of the domain X is represented as a N-variable state vector.If one waits long enough after applying small external stimulus h at time t ¼ 0, the only persisting effect will be a net shift in the spiral's collective coordinates, i.e.
for some linear functional W l .When working in the continuum approximation of the medium (i.e., X & R 2 , the linear functional acts as an inner product with hfjgi ¼ Ð R 2 f H gdxdy and : H the Hermitian transpose.The functions W l are known as the "response functions" (RFs) of the spiral wave solution. 12In essence, every point of a response function for a given state variable records the spiral wave's drift response to a perfectly localized stimulus applied to the same state variable.Using h j to denote the j-th state variable component of h, we have In that sense, spiral wave RFs are similar to the impulse response in engineering or Green functions in physics.Moreover, the phase component W W is the spatial generalisation of the phase response curve (PRC) for oscillating systems.Knowledge of the spiral wave response functions enables predicting the spiral wave drift in the perturbative regime.Overlap integrals of spiral wave RFs thus appear in the equations of motion for three-dimensional scroll wave filaments, [13][14][15][16] or in the theory of 2D spiral waves drifting due to a constant external field, 17 surface curvature, 18 or mechano-electrical feedback. 19In one spatial dimension, the translational RF of the wave front determines the velocitycurvature relation 20,21 and the shape of the RF itself can be used to shape reaction-diffusion patterns. 22f the underlying model equations are known (e.g., for a cardiac tissue model), the response functions for a (meandering) spiral can be found by linearising the system around its relative equilibrium (periodic orbit).The critical adjoint eigenfunctions to the associated linear operator are precisely the response functions; see e.g., Refs.23-25 for details.This way, spiral wave RFs have been numerically computed, first for simple phenomenological models with few state variables 12,26,27 and later with increasing accuracy, 28 up to the Beeler-Reuter cardiac tissue model. 29In all cases, the spiral wave response functions turned out to be exponentially localised in the core region, and in the view of Eq. (3), this property grants spiral waves their so-called "particle-wave dualism." 12 Only for the complex Ginzburg-Landau equation, some asymptotic expressions were found. 26or meandering spiral waves, the first RFs have only been computed very recently by Marcotte and Grigoriev 30 by solving the adjoint linear system.In the case of meander, the spiral wave solution and its response functions are periodic in the meander phase W in a slowly rotating frame.A challenge that was surmounted in the numerical computation of Ref. 30 was that the full spiral solution (with period 300 ms) and typical discretisation time step (0.1 ms) do not fit into 8 GB of RAM memory and the solution needs to be iterated many times from the initial condition ðW ¼ 0Þ over the meander period to compute the response function.For a 3-variable atrial model, the RFs as well as the first tens of leading eigenfunctions were computed in 72 h, using graphical processing units to accelerate the time stepping in the forward problem.
In this work, we present an alternative approach to evaluate spiral wave RFs.The idea is to use Eq. ( 3) directly and perturb an initial condition in different simulations with a set of elementary responses and record the resulting shifts in the collective coordinates.
Section II describes the numerical methods and models used in the forward evolution, and subsequent RF reconstruction.In Sec.III, we present the RFs for a circular-core and meandering spiral in Barkley's model, and a linear core spiral in the Fenton-Karma (FK) model.

A. Theory
We consider spiral waves as particular solutions of a non-linear dynamical system in 2 spatial dimensions x, y with N state variables _ uðx; y; tÞ ¼ Q uðx; y; tÞ ½ þhðx; y; tÞ; where uðx; y; tÞ : The operator Q is taken to be space and time invariant (disregarding boundary conditions at a large distance).A common choice to model oscillatory and excitable media is to take a reaction-diffusion system, i.e.
A spiral wave solution solves (4) in a suitable comoving frame of reference, with drift speeds V x ðW; UÞ; V y ðW; UÞ, rotation speed x 0 V h ðWÞ, while its meander phase varies as _ W ¼ V W X 0 .The unperturbed spiral wave solution u 0 ðx; y; t; X l ðtÞÞ thus satisfies The second equation describes the unperturbed spiral tip trajectory.
Ideally, one could subject the system to a stimulus (3) that is localised at a single grid point, but for small stimulus amplitude , the resulting shifts DX l will be below the measurement threshold in the discretised system.Therefore, we take a broader stimulus whose profile is a Gaussian kernel function, in the k-th state variable We choose to be equal to 0:02 À 0:05 times the wave amplitude, and then choose r large enough to make the coordinate shifts measurable in the system.Typically, r was between 2 and 10 times the simulation grid resolution, see Table I.We are also free to choose the locations of where the external stimuli are centered, and took them on a rectangular grid of spacing r and size N x;stim Â N y;stim around the position of the spiral wave tip at the time of stimulus delivery.The instantaneous delivery of a stimulus, i.e., dðtÞ was implemented by adding the stimulus at one time step in the simulation.As a result, when the perturbation amplitude is large enough to dominate discretisation effects of the lattice, the shifts D l will be the spatial convolution of the response function with the Gaussian kernel Here, k is the index to the state variable for which the RF is currently probed.However, if the kernel G (in simulations) or an estimate to it (in experiments) is known, one can perform a deconvolution to the measured shifts in order to reconstruct the RF.
To allow the possibility of an unstructured grid of measurement points, we deconvolve using Tikhonov regularisation 31 rather than a Fourier-based method.Subject the system in M experiments to M perturbations centered at the positions ðx i 0 ; y i 0 Þ and record the shifts b i ¼ ðDX l Þ i =.We zero-padded the grid of responses b i with 5 rows in x and y directions.
When reconstructing the RFs in N points ðx j r ; y j r Þ, this leads to a M Â N linear system Here, we choose the reconstruction grid equal to the stimulation grid.We make A sparse by rounding an element to zero if its value lies below 0.001.As a regularizing term, we take the 5-point Laplacian stencil of the set of reconstruction points, denoted by a matrix C. Thereafter, we find the response function using We manually choose the lowest possible k which still suppresses noise outside the spiral core region, see values in Table I.The inversion is performed by the mldivide command in matlab 32 using Cholesky factorization.

B. Reaction-diffusion models
For definiteness, we illustrate the method for spiral waves in a reaction-diffusion finite-difference model.It can also be applied in other systems supporting spiral waves (or, more general, relative periodic orbits), as long as Euclidean symmetry of the system is well enough preserved.Factors that may break this symmetry are the underlying discrete nature of the medium (e.g., biological cells or a computational grid), spatially varying anisotropy, and inhomogeneities in an experimental preparation or tissue.
In this work, we illustrate our methods using numerical simulations only.First, we use a reaction-diffusion model with Barkley's kinetics, 33 i.e., u ¼ ½u v T , P ¼ diagð1; 0Þ and We will take the following parameters sets: which, respectively, produce a circular-core spiral wave and a meander spiral wave with flower-like tip trajectory, see Fig. 1.We found the tip position every 0.1 time units as the intersection of the isolines u ¼ 0.5 and v ¼ 0:5a À b.
Secondly, we will determine the RFs for a cardiac tissue model with a star-like tip trajectory, which is sometimes called a "linear core."Hereto, we take the three variable Fenton-Karma (FK) model with guinea pig (GP) parameters. 11It has u ¼ ½u v w T and P ¼ diagð0:1; 0; 0Þmm 2 =ms.While u is the normalized transmembrane potential, v and w are fast and slow recovery variables.The details of the phenomenological electrical currents in the model can be found in Ref. 11.The model produces a spiral wave with meander period 75 ms, as shown in Fig. 1(e).We determined the tip position as the intersection of the isoline u ¼ 0.5 with the same curve 1 ms before.
The reaction-diffusion models were numerically integrated with the explicit Euler scheme one a Cartesian grid with Neumann boundary conditions, see Table I.

C. Measurement of the shifts in collective coordinates
From Eq. (1), the response function values can be inferred from shifts in the collective coordinates, i.e., the center, rotation angle, and time phase of the meander tip trajectories.It is thus necessary to determine accurately those shifts, at a resolution that is finer than the lattice size in simulations.
Let us first treat the circular core case, see Fig. 2(a).After applying a perturbation to a spiral wave solution at time t ¼ 0, we let the system evolve for several periods, and during a time interval ½t start ; t end , we sample the tip trajectory every dt tip .The i-th measurement is taken at time t i 2 ½t start ; t end , yielding corresponding tip positions (x i , y i ).These points will lie approximately on a circle (see Fig. 2 left), whose center (X, Y) and radius r can be found by linear regression. 34Thereafter, the instantaneous rotation phase U i can be found as the polar angle corresponding to (x i , y i ).
Finally, we perform a least squares linear regression The case of meander requires more fine-tuning to the given tip trajectory, as shown in Figs.2(c) and 2(d).For the BA2 model with the chosen, we first list all self-intersections of the tip trajectory.From this list, we remove self-intersections for which the time between the first and second visit is longer than 10 time units.We exclude longer time intervals, since the tip trajectory shown in Fig. 2(b) will reach again the same region after the completion of the meander flower.Such intersection points are not useful to determine the meander phase.
Since the motion is quasi-periodic, the time difference between subsequent intersection points will repeat itself after N rep points, with N rep an integer that is small if one is far from the onset of meander.To find N rep , we looped over the possible number of classes N rep between 1 and 10, and computed the standard deviation at times between occurrence of a point, and found that it nearly vanished for N rep ¼ 6 classes of intersection points.In Fig. 2(b), it can be seen that there are 3 classes of intersection points when looking at the trajectory as a whole, but they are all visited twice, and we find N rep ¼ 6.As each class of points lies on a circle, we estimate the centre and radius of those circles and then choose the class with smallest radius to be the set of fiducial points to determine the collective coordinates, see Figs. 1(b) and 2(b).On this circle of fiducial points, we repeat the procedure for circular-core spiral waves to find X; Y; U 0 and x 0 .At the j-th fiducial point, we also take W j ¼ j2p and then fit Thirdly, we consider the linear-core example with FK kinetics, as shown in Fig. 2(c).As the fiducial points, we will take the turning points of the trajectory, which also lie on a circle if the meander pattern is quasi-periodic.First, we smooth the tip trajectory between t 1 and t 2 using a Gaussian window of size 25.Then, we take the average of the lists x i , y i as the first estimate of X, Y. Next, we determine the points where the distance to (X, Y) reaches a maximum and fit a circle through them; the centre of this circle is an update estimate of (X,Y).After repeating this iterative process twice, (X, Y) remained constant and were taken as the meander centre position.Thereafter, the turning points were taken as fiducial points and the procedure for flower-like meander described above was used to determine U 0 ; W 0 ; x 0 , and X 0 .

A. Circular-core spiral with Barkley kinetics
For the BA1 model, we computed the three RFs using (12) and compare them to the results of solving the adjoint method in Fig. 3.For the latter, we used the publicly available software dxspiral, 28,35 which was developed to compute the RFs for reaction-diffusion models that support rigidly rotating spiral waves.From Fig. 3, it can be noted that the final result using our method closely matches the result obtained using dxspiral in amplitude and shape of the profiles.For the RF measurement, we sampled the Gaussian center every two grid points of the computation grid, i.e., with spatial resolution 0.2, to obtain the RF in a set of 101 Â 101 points.As every forward run of 80 time units took 22 s on our system (cþþ program using MPI domain splitting on 16 CPUs), the total computation time amounted to 10.3 h per variable, which is two orders of magnitude slower For the circular-core case [panel (a)], points on the tip trajectory itself lie on a circle from which the rotation phase can be regressed.For the meander cases (b) and (c), a circle is fitted to the fiducial points (red) in the tip trajectory in order to find the position of the meander centre (X, Y) (leftmost panels, indicated by þ).The polar angles U i of the fiducial points with respect to (X, Y) are regressed vs. time (top right panels) to find U 0 and x 0 .For meandering spiral waves (BA2, FK), the times at which the fiducial points are visited are used to compute the absolute meander phase W 0 and meander period T 0 ¼ 2p=X 0 (bottom right panels).
than the time of 5-10 min.that is needed using dxspiral.Nevertheless, by this comparison we show that the measurement method can be used to compute RFs in detail.The difference between both methods is displayed in Fig. 4.

B. Flower-like meander with Barkley kinetics
For the BA2 model parameters, the obtained RFs are shown in Figs. 5 and 6.Note that in a frame of reference with rotation frequency x 0 , the RFs would be periodic in W.
FIG. 3. Comparison of computed RFs using the adjoint method (columns 1-2) and measurement method (columns 3-4) for a counter-clockwise circular-core spiral wave with Barkley kinetics BA1.The black line is the isoline u ¼ 0.5; the spiral wave tip is situated in the top-left quadrant.Results in columns 1-2 were computed with dxspiral on a polar grid N r ¼ 120, R ¼ 12, N h ¼ 128.Results in columns 3-4 were found using the measurement method with dx ¼ 0.1, N x ¼ N y ¼ 240, dt ¼ 0.002375, tip trajectory fitting for t 2 ½40; 80, tip calculation every 0.1 t.u., A ¼ 0.05, k ¼ 1=ð2dxÞ, k ¼ 10.FIG. 6. Same as Fig. 5, at a different phase of the meander cycle, i.e., 2.5 t.u. or 0.5 meander periods later.
We show two cases of W in Figs. 5 and 6.Here too, the four RFs are strongly localized in the tip region.
C. Linear-core meander with Fenton-Karma kinetics Applying the measurement procedure to the FK-GP model delivers the RFs shown in Figs.7 and 8.For these two frames in the meander cycle, one can see that for all state variables, the RFs are localised near the turning points of the tip trajectory.If this holds for all times during the meander cycle, it is an important finding, since it means that external stimuli cannot change the rotor's path if they are not applied at locations where the few next turning points will occur.The FK model was created to represent fast and slow recovery processes by the v and w variables, respectively.From the RF structure, we note that the W l v are solely peaked at the next turning point, which can be understood since a stimulus in v will have decayed when the second turning point is reached after 2T 0 ¼ 150 ms.Contrarily, perturbing w will affect significantly the spiral wave dynamics at the 4 upcoming points, as is evident from the W W w and W U w components.As expected, the effect on the phases W; U has the same sign for all the turning points.This is not the case for the W x and W y components, since they depend on the orientation (i.e., have a tensor index which needs to be transformed).There, we also plot ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi for all 3 state variables in Fig. 9.For the rotational modes too, we note the  memory effect is longest for the w component.Moreover, the maximal amplitude is only found at the second turning point, probably since the tissue there is unexcited at the time of stimulus, contrary to the location of the first turning point.

IV. DISCUSSION
We have developed and presented an "experimental" method to probe spiral wave response functions, and validated it on a numerical simulation of a circular-core spiral wave.In terms of numerical simulations, the method used is rather inefficient compared with other existing methods, but we nevertheless want to demonstrate three important points.
First, if one has sufficient computational power at one's disposal, it is possible to compute spiral wave RFs using minimal adaptation to the existing computer codes.By performing runs sequentially or on a cluster, there are no bigger memory limitations than for solving the original system.
Second, we see most use of the presented method in a truly experimental determination of spiral wave response functions, e.g., in the BZ-reaction or cardiac monolayers where Euclidean symmetry is more likely to be observed than in cardiac tissue.Note that in real-life experiments, the preparation only needs to be set up once.For, if the medium properties are sufficiently homogeneous (e.g., in BZ-reaction), one can perturb the spiral, let it rotate for 5-10 rotations to measure the shifts, then perturb again and so on.As a result, one may, e.g., infer 20 points of the RF from 100 spiral wave rotations.For cardiac tissue, pinning effects may be pronounced and therefore one may turn to monolayers of cultured tissue first.Some future applications may lie outside the context of reaction-diffusion systems.The location and timing of the stimuli need not occur on a regular grid in time or space, enabling to take a Monte-Carlo approach, and better determine the RF at every new experiment in the series.
Third, by computing the RF for the FK-GP cardiac tissue model, we have demonstrated that the spiral wave RFs are concentrated near the turning points of the tip trajectory.Thus, the efforts to control or steer rotors in cardiac tissue, even if delivered to the entire tissue, are only effective near the turning points.Conversely, the interaction of a cardiac rotor with ambient structure (inhomogeneities, anisotropy, and tissue thickness) in quasi-2D and 3D tissue is also likely to be dominated by the structures it encounters at the position of its turning points.
FIG. 1. Spiral waves and tip trajectories for the three reaction-diffusion models used, in Barkley's model with circular core [BA1, panel (a)] and meandering core [BA2, panel (c)] and the FK-GP model [FK, panel (e)].Panels (b), (d), and (f) show how the collective coordinates X; Y; U; W were fitted using fiducial points (red dots) within a meander or rotation cycle.

FIG. 2 .
FIG. 2. Measurement of collective coordinates in the BA1 model (a), BA2 model (b), and FK model (c).For the circular-core case [panel (a)], points on the tip trajectory itself lie on a circle from which the rotation phase can be regressed.For the meander cases (b) and (c), a circle is fitted to the fiducial points (red) in the tip trajectory in order to find the position of the meander centre (X, Y) (leftmost panels, indicated by þ).The polar angles U i of the fiducial points with respect to (X, Y) are regressed vs. time (top right panels) to find U 0 and x 0 .For meandering spiral waves (BA2, FK), the times at which the fiducial points are visited are used to compute the absolute meander phase W 0 and meander period T 0 ¼ 2p=X 0 (bottom right panels).

FIG. 4 .
FIG. 4. Difference between the computed RFs using the adjoint method (dxspiral) and measurement method for the spiral wave from Fig.3.

FIG. 7 .
FIG. 7. Computed response functions for the FK cardiac tissue model, for a given meander phase denoted by spiral wave position (black line, u ¼ 0.5) and tip position (red dot).The white trace is the upcoming unperturbed tip trajectory.

TABLE I .
Parameters used in numerical methods for the different models.