A theory for spiral wave drift in reaction-diffusion-mechanics systems

Reaction-diffusion mechanics (RDM) systems describe a wide range of practically important phenomena where deformation substantially affects wave and vortex dynamics. Here, we develop the first theory to describe the dynamics of rotating spiral waves in RDM systems, combining response function theory with a mechanical Green’s function. This theory explains the mechanically-induced drift of spiral waves as a resonance phenomenon, and it can predict the drift trajectories and the final attractors from measurable characteristics of the system. Theoretical predictions are confirmed by numerical simulations. The results can be applied to cardiac tissue, where the drift of spiral waves is an important factor in determining different types of cardiac arrhythmias.


Introduction
Spiral waves have been observed in a wide variety of physical, biological, and chemical systems far from equilibrium. Examples include chemical waves in Belousov-Zhabotinsky (BZ) reactions [1], waves of CO oxidation on platinum surfaces [2], cyclic adenosine monophosphate (cAMP) waves during amoeba morphogenesis [3], waves of spreading depression in the brain [4], and electrical waves in cardiac tissue [5,6]. Furthermore, spiral waves have exhibited very rich and complicated dynamics due to the intrinsic properties of the media or in response to external fields [7,8]. The drift of spiral waves is among the most prominent examples that has been observed in BZ chemical systems [9] or studied in biological systems such as the heart [10,11]. In the heart, the drift of spiral waves is believed to be responsible for certain life-threatening cardiac arrhythmias [12].
The propagation of spiral waves is typically accompanied by other important processes such as the mechanical deformation of the medium [3,13,14]. In the heart, propagating electrical waves initiate cardiac contraction, which in turn affects their propagation, a process known as mechanoelectrical feedback (MEF) [15]. Although deformation is known to be crucial in the mentioned systems, most previous theoretical and experimental studies have not addressed the combined effects of the medium mechanics and spiral wave dynamics. In cardiac tissue, the most immediate, and perhaps most important, effect of this MEF is the activation of stretch-activated depolarizing currents [16], which may underlie the initiation of cardiac arrhythmias [15].
To model the basic MEF effects on wave propagation, the concept of coupled reaction-diffusion mechanics (RDM) systems [17][18][19] was proposed; this concept combines a very general description of deformation with classical reaction-diffusion kinetics. A series of papers has shown that this kind of MEF induces complex dynamics of waves, such as drifting pacemakers [18,20] and the break-up [19], drift [19,20], initiation [20,21], and unpinning [22] of spiral waves. However, in those cases the results were obtained using numerical simulations for particular reaction kinetics and mechanical models only. No theories have yet been developed to describe the dynamics of spiral waves in RDM systems using analytical approaches. Such analytical theories are believed to be important tools that can help us find out how general the numerical findings are; they can also pave a way for applications to practically important problems.
In this work we present an analytical approach to study spiral wave dynamics in RDM systems. We combine response function theory [23] with a Green's function formalism to account for mechanical influence, and we derive an equation for spiral wave drift. The MEF-induced drift is reduced to a resonant forcing problem in the rotating frame of the spiral: during its rotation, the spiral wave perceives a time-varying perturbation since the domain boundaries are not stationary in the spiral's frame of reference. The resonant component of this boundary-induced forcing yields the net spiral drift.
Our theoretical predictions are compared to numerical simulations. We find the relative angle and magnitude of the spiral wave's drift and identify the spatial attractors of the system. Some of them, such as the center of the medium, have already been reported in numerical simulations [19]. Using this theory, we also find new regimes and attractors that also are confirmed by simulations. For example, our theory predicts that the center of the domain can also be repulsive, and that multiple stable dynamical attractors may coexist. Although different from previous findings, these predictions are confirmed by numerical simulations.
The developed analytical approach allows us to generalize the numerical results, as our analytical findings are based on an Archimedean description of spiral wave geometry, which is common for all types of excitable media.

Model
The reaction and MEF parts of our model for cardiac tissue are as in [18,20,24,25], supplemented here with the Navier-Cauchy equilibrium equations from linear elasticity to facilitate analytical calculations. In a twodimensional medium with Cartesian material coordinates (x, y), the transmembrane voltage u and recovery variable v are evolved according to modified Aliev-Panfilov kinetics [26], as used in [18]: We set ϵ = u ( ) 1.0 when < u a, and 0.1 otherwise. Unlike some previous works [18,25], the discontinuity of ϵ u ( )occurs here at u = a, not at a fixed value of u = 0.05. In numerical simulations, we take = = k k 8, 1.5 T , and the model parameter a will be varied to change wavelength and excitability. Active mechanical tension T is developed in our model as [17,19,24,25] T We model the tissue as isotropic and elastic, and we furthermore assume that spatial displacements are small. Then, the Lagrangian displacement, where ⃗ x lab denotes Eulerian (laboratory) coordinates and ⃗ x are Lagrangian (material) coordinates. We assume that the displacement ⃗ U instantaneously obeys the Navier-Cauchy equations for mechanical equilibrium: with Lamé parameters λ and μ. Elastostatic conditions have also been assumed in previous electromechanical studies of MEF in the heart [17,18,24,25,27]. In previous work using large displacements [19], the stretch, S, was equal to − C 1, where C is the determinant of the right Cauchy-Green deformation tensor. Presently, we work in the regime of small deformations, where In our model, the feedback of deformation on excitation is given by stretch-activated membrane currents, I , s in (1). Experimental studies have shown that stretch-activated channels immediately respond to mechanical stretch and follow a linear current-voltage relationship [16,28]. Linear models for I s have been proposed [29,30]. Therefore, we choose the linear current-voltage relation for the stretch-activated current in (1), as used in [18,24]: Equations (1)-(6) form a closed feedback loop for mechanoelectrical coupling in the domain, Ω. As is customary, we impose Neumann boundary conditions on the state variable u at the edge, Ω ∂ , of the domain. For the mechanical subsystem, we work with fixated boundaries-that is, The same boundary condition was applied before in [18,19]. We use this assumption to mimic the isovolumic phases of the cardiac cycle and isometric boundary conditions in tissue experiments.

Material coordinates
We have taken a Lagrangian standpoint here: All spatial dependencies in (1)-(6) refer to two-dimensional material coordinates, ⃗ x , which we will also write as x a ∈ a ( [1, 2])in index notation. In material coordinates x a , we use = to denote the position of the spiral's rotation center at each instance of time.
Laboratory coordinates ⃗ x lab will only be used in section 6.1 to quantify the physical displacement, ⃗ − ⃗ x x lab , due to MEF.

Working assumptions
In our theory, we will linearize the equations around both small deformations and the rotating spiral wave solution. Additionally, the sensitivity of spiral waves to external stimuli is assumed to be strongly localized, such that response function theory [31][32][33] can be applied. Thus, we make the following three assumptions: , with η being small). We pursue leading-order dynamics in η.
(ii) In the absence of mechanical feedback ( = I 0 s ), (1) and (2) can be written as an evolution for the state = u v u [ ] T . We assume there exists a rigidly rotating spiral wave solution, = u v u [ ] T 0 0 0 , to the RD system, (1) and (2), with rotation frequency ω 0 . The solution u 0 is stable with respect to small perturbations (i.e., it is an attractor in phase space). We will denote the unperturbed spiral solution with rotation center, (X, Y) and rotation phase ϕ as ϕ ⃗ ⃗ x X u ( ; , ) 0 . In terms of polar coordinates around the rotation center, sin , the rigid rotation condition can be expressed as When a spiral wave is subjected to a small spatiotemporal perturbation, ⃗ x t h( , ), linear superposition can be applied due to (i), and by (ii) the system quickly relaxes towards u 0 for some position (X, Y) and rotation phase ϕ. Thus, in an unbounded domain, the net spiral motion at a given time is a spatial convolution over all sources: Here, the functions a a a u v and a sum over state variables u v , is implied in (8). The functions W a , θ W are known as translational and rotational response functions [33] and have been used in many quantitative descriptions of spiral and scroll wave dynamics [11,31,[34][35][36][37][38][39][40]. Our third hypothesis is that W a and θ W exponentially decay with the distance from the spiral core region, with spatial constant d core , which is indicated in figure 1. This property was proven for a specific model [33] and has been verified numerically for reactiondiffusion models with few state variables and smooth reaction kinetics [36,[41][42][43]. However, general conditions for the localization of response functions currently remain unknown.
Note that the noncontinuously differentiable reaction term in (2) does not forbid the use of response functions. However, this discontinuity will yield a delta-distribution term in the Jacobian in places where the associated response functions will be discontinuous but nonsingular. Since our results will involve overlap integrals using response functions, those integrals are well-posed and our theoretical results are therefore also valid for noncontinuously differentiable reaction kinetics. However, the numerical evaluation of response functions in this case is cumbersome due to the singular Jacobian. Therefore, we will not attempt to evaluate the response function for the reaction-diffusion model (1) and (2), but we will compare theoretical predictions of relative drift magnitude and direction with the outcome of numerical simulations.
Hypothesis (iii) is an assumption on the electrical subsystem (1) and (2). In the absence of mechanical deformation, it predicts that spiral-boundary interaction causes drift of the spiral that attenuates exponentially with the distance, d, to the boundary (i.e. ∂ ∝ − X dd exp( ) t a core ). However, we will show later in this paper that a mechanical boundary with a no-displacement condition at a distance, d, induces a stretch-activated current proportional to d 1 2 . Therefore, if the spiral core lies further than d core from the boundary, the electrical perturbations of the pattern cause an exponentially small drift velocity, while the MEF induces drift of order 1. We will work in this regime, and accept that our theory will not be valid if the tip of the spiral comes closer than d core to the domain boundary.
If the spiral wave's rotation center is located at a distance, > d d core , from the boundary, the spiral profile will still closely resemble ϕ ⃗ ⃗ x X u ( ; , ) 0 inside the domain. However, near the boundary, the Neumann condition on u imposes that the wave front should intersect the domain boundary orthogonally. Therefore, a truncation error is induced in the wave profile in a layer of thickness, d trunc , close to the boundary. From the velocity-curvature relation for wavefronts [44,45], it follows that d trunc is of the order of the wave-front thickness, and also of the same order as the critical radius for excitation. We will work in the regime where In what follows, we will ignore the domain truncation error.

Spiral wave drift equations
Provided that the unperturbed spiral wave solution is a phase-space attractor by condition (ii), a dynamical spiral wave state will, apart from its position ⃗ X and rotation phase ϕ, quickly lose its dependence on previous states. Therefore, in the domain Ω, the spiral wave evolves as For the mechanically coupled system, we use response function expressions (8) to find the resulting spiral wave drift. From (1), we have u s Note that the response functions, μ W u , would be constant in a frame attached to the spiral rotation center that rotates at the spiral frequency. Therefore, they depend on time only via ⃗ X t ( )and ϕ t ( ). The same dependency will be shown to hold for the stretch-activated current: Given an unperturbed spiral wave, Since the Navier-Cauchy equation (4) are linear, the resulting stretch can be found using a mechanical Green's function: U satisfying (4) and the boundary conditions (7) is denoted as Taking the divergence of both sides and integrating by parts delivers , given in (6), we have established that I s only depends on time through ⃗ X t ( )and ϕ t ( ). Returning to (10), we have justified the dependency of time in the right-hand sides of (3.3).

Average spiral drift
Under low MEF, the (3.3) can be averaged over one rotation period. When the spiral is far away from the boundaries ( ≫ d d core ), only the MEF induces spiral drift, which is therefore of order η. Within a rotation period of O(1), the spiral's rotation center thus has Then, (9b) becomes an ordinary differential equation, which can be integrated to Here the integral of g(t) vanishes over one period, and The net spatial spiral drift over one period can be found by: Since the f a are π 2 -periodic functions of ϕ, they can be expanded in a Fourier series a n n a n i Within one rotation period, we can take a a n n a n n a n Thus, within one rotation period, the harmonics in ϕ ± e n i describes an epicycle trajectory for the spiral wave's rotation center. However, only C a 0 contributes to the net spiral drift over one rotation period. If one averages spiral position over one period, its mean velocity is a a 0 0 2 Substituting the explicit expressions (12) and (10), we arrive at following expression for the period-averaged drift of spiral waves under MEF:

Net spiral wave drift in the corotating frame
In the previous paragraph, we found that net spiral drift is given by the vector ϕ ⃗ ⃗ f X ( , )averaged over all rotation phases. Looking at the right-hand side of (19), we may now treat the spiral position ⃗ X as a formal parameter, without explicit time dependency.
Until now, we have used material coordinates ⃗ x in a nonrotating frame, which we denote as x a , ∈ a {1, 2}. With its rotation center at a given location of the medium, a spiral wave rotates at an average frequency, , 2}) in index notation. These are related to the non-rotating material coordinates, x a , as Since a spiral wave's sensitivity to external stimuli is time-independent in the corotating frame, in the corotating coordinates x A , the response functions are time-independent functions Due to the localization of response functions, the integration domain can be limited to the spiral core region, Ω c . The net spiral drift becomes Equation (23) can be rewritten in complex notation: Defining complex velocity = + x y , and a complex-valued response function, This expression shows that the net MEF-induced spiral drift results from the first Fourier component of the induced stretch distribution. This Fourier component is multiplied with a weighting function, , evaluated at S = 0, and integrated over the tip region to produce the net spiral drift. In terms of active tension development, combining (19) and (24) delivers Furthermore, we note that the mechanical Green's function slowly varies on the scale of the domain size, L 2 , which is much larger than the decay width, ξ, of the translational response function, + W u . Thus we may write in the lowest order in ξ L that Higher-order corrections will contain partial derivatives of G S , evaluated at the spiral tip. Hence, (25) can be written as the product Expressions (26)- (28) show that the magnitude, = | | + V V , and direction, + V arg( ), can be computed from the first Fourier component of the induced stretch at the spiral tip, up to a constant amplitude, H, and absolute drift angle α. In section 6 we will compare (27) with the drift resulting from numerical simulations.
3.6. Time-course of active tension development If the radius described by the spiral tip is not large compared to the domain size, L 2 , and the spiral is far from the domain boundary, the induced stretch does not vary much in the spiral's core region for a given rotation phase. In this case, we can approximate ; .
We now show that the resulting spiral wave drift is largely independent of the detailed time course of active tension development in the tissue; only the delay between excitation and active tension development will play a role. In the electrical subsystem (1)-(3) and in the limit of small MEF, the active tension develops regardless of the local stretch and its derivatives. Instead, the active tension only depends on the time elapsed since the point was electrically excited, or, equivalently, its excitation phase, Φ. At any given time, a counterclockwise rotating spiral wave with wavefront shape θ Θ = r ( ) 0 in polar coordinates θ r ( , ) around the its rotation center has In material coordinates Expressing that the active tension developed in our model only depends on the excitation phase yields, for points outside the spiral core region: 0 a c t 0 Here, the function Φ T ( ) act is the temporal profile of the periodically developed active tension in each point of the unperturbed spiral wave, u 0 , outside its core region. Furthermore, we show below in section 4 that if the spiral core is not close to the domain boundary, maximal stretch occurs near the boundaries, which is thus outside the spiral core. Then, after applying (34), we may change integration variable Here, a new shape factor, A T , appears, which equals the first Fourier coefficient of the time-course of active tension development: T i˜a ct Expression (35) shows that, within the approximations made, the stationary positions for spiral waves in the RDM system (1)-(6) do not depend on the details of the development of active tension over time in each active element (cell) of the medium. Changing the function Φ T ( ) act only changes the complex prefactor A T in (35). However, the time-course of active tension development may alter the stability of the equilibria, and its effect can be computed by computing A T . If, for example, active tension is delayed over a time t 1 , a complex pre-factor ω e t i 1 is added to + V , which can change dynamical attractors to repulsive sites, and vice versa. In conclusion, (35) reveals that of the electrical subsystem, only the shape of the spiral wave's front determines the position of the attractors. This property can be made explicit by considering a reduced system in which tension is instantaneously developed at the wavefront (i.e., with πδ ϕ = T 2 (˜) act ). By (36), the drift velocity, + Ṽ , in this system is equal to ⃗ + V X ( )of the original system, divided by A T . However, putting in the instantaneously developed tension already in (30) would reduce the surface integral over the domain to a line integral. Equating both expressions for + Ṽ in the reduced system allows one to rewrite the drift velocity in the original system in terms of a line integral over the wavefront surface: is a parameterization of the wavefront surface, , given its rotation phase, ϕ. Given that the spiral waves that form in reaction-diffusion systems asymptotically tend to an Archimedean spiral shape, the spiral's wavelength, Λ, is the only characteristic of the electrical subsystem that determines the position of attractors in the RDM system.

Position and stability of equilibria
We are interested in identifying the equilibria and their linear stability given the complex-valued velocity field, . Linear stability analysis of a two-dimensional system shows that the simple (i.e., linear) equilibria can be classified into stable nodes, stars or spirals, unstable nodes, stars or spirals, saddle points, and center points [46]. Example trajectories and the complex velocity field for six cases are shown in figure 2. In this work, we will not discriminate between nodes, stars, and spirals. Stable nodes, stars, and spirals will be denoted as attractive equilibria, while we refer to unstable nodes, stars, and spirals as repulsive equilibria.
We note that all types of linear equilibria correspond to phase singularities of the velocity field, + V X Y ( , ). To locate phase singularities numerically on a Cartesian grid, we track the complex phase, ψ = + V arg( ), counterclockwise around every cell of the grid, and see whether it changes by π ±2 . If so, a phase singularity is present, and we locally fit a linear system, 22 . Then, it is well known that the trace and determinant of the coefficient matrix A determine the type of equilibrium [46]:

Mechanical Green's function for half-plane and square domains
We will compare our general theory above with numerical simulations on a square domain of size = L 2 60. To make predictions on spiral wave drift, we need to know the mechanical Green's function, ⃗ ⃗ G x x ( ; ) S 0 , for the Navier-Cauchy equations in a given geometry, or an approximation of it.
The Green's function for the half-plane > y 0 with a zero-displacement boundary condition at y = 0 can be found analytically. Using the method of complex variables in elasticity, as outlined in [47], we find the following fundamental solution to the Navier-Cauchy equations for stretch  = ⃗ ⃗ S U · induced by active tension,    An example of the resulting stretch distribution is displayed in figure 3(b).

Electromechanical simulations
We performed numerical simulations using the discrete reaction-diffusion mechanics (dRDM) model developed in [24]. The dRDM model constitutes a generic electromechanical model for cardiac tissue using (1) and (2) to describe the excitation processes, and (3) to model the excitation-contraction coupling. The model is based on representing tissue by mass points connected by active and passive springs. The MEF of cardiac tissue is modeled in dRDM by stretch-activated currents, I s , with an equation similar to (6), which describes stretch not in terms of the deformation tensor, but of the surface area of a quadrilateral formed by neighboring mass points [24]. We used = = = = k k g E 8, 1.5, 1, and 1 T s s as parameters for solving the modelʼs equations. Furthermore, the model uses a mass-lattice system whose constitutive relations can be approximated with the Seth material relation [24].
We study a model with a square domain of size Thereafter, we averaged tip positions over one spiral period. Hereto, we deduced rotation phase ϕ as the angle between the tangent vector to the tip trajectory and the positive x-axis, and we added multiples of π 2 to obtain a monotonous function. To find the average tip position at a time with phase ϕ 0 , we take the barycenter of the segment of the trajectory that has phases between ϕ π ϕ π − + and 0 0 . For given model parameter a, we interpolated the Cartesian components, V V , , x y of the period-averaged drift velocities onto a finer rectangular grid, and we extracted the complex velocity phase, ψ = V V arctan( ) y x . The results are displayed in the top row of figure 6.

Theoretical prediction of drift velocity
For comparison, we also used the theoretical drift expression (30) and the mechanical Green's function (39) to predict the relative magnitude and angle of spiral wave drift + V rel using (30). Hereto, we first characterized the spiral wave by measuring its wavelength Λ, core radius, and periodically developed active tension,  (39) to yield stretch at the tip position, after which negative stretch values were set to zero to include the ramp function, P, from (6). Finally, the first Fourier coefficient with respect to the rotation phase was computed to yield the complex-valued relative drift velocity, + V rel . From (30), this quantity is expected to relate to the absolute drift velocity, = i rel for real-valued constants, α H, . The constants α H, that determine the absolute magnitude and direction of drift can be evaluated using (28) if the response functions are known. Existing tools to compute spiral wave response functions (e.g., [42]) require continuously differentiable reaction kinetics. In the current study, however, noncontinuously differentiable reaction kinetics were used to maintain compatibility with previous works. Therefore, we fitted the absolute drift angle, α, in (5.2) from the observed drift angle at the start of each dRDM simulation. We did not fit the magnitude factor, H , since it affects neither the stability of equilibria, nor the spiral wave's drift trajectory.
With α fitted from initial drift velocity, we obtained from (5.2) a velocity field, ⃗ V X Y ( , ), which was numerically integrated using Euler stepping to produce the theoretical drift trajectories shown in figures 5 and 6.

Drift trajectories
The results of the dRDM simulations and their comparison with theoretical predictions are shown in figure 6. Due to rotation symmetry over π + n 2, 4 1 equilibria are found in each case with n integer. The results shown are for a clockwise rotating spiral wave. For the opposite chirality, drift trajectories will be mirrored.
The results from theory and numerical simulations are qualitatively similar in terms of the number and position of attractors and the shape of the spiral drift trajectories.
For large = a 0.08 (i.e., long spiral wavelength and low excitability), a single attractor of the spiral type is found. When a is decreased to = a 0.07, both theory and simulations show four new attractors close to the corners of the domain. Closer investigation of this transition shows that close to = a 0.076, the central attracting spiral becomes unstable through a Hopf bifurcation, yielding nearly circular limit cycles, as shown in figures 7(d)-(f).
When a is further reduced to = a 0.066, another Hopf bifurcation occurs; see figure 7(b). At this critical value, the four corner attractors lose stability while the square's midpoint stabilizes, as seen in figure 6 for = a 0.06. No qualitative changes take place when a is further lowered to 0.05.
While the absolute drift direction (i.e., the angle α) is fitted between the theoretical and simulation panels of figure 6, the position of equilibria results from the integral in (30) alone. From symmetry, the square's center is always stationary, which is also seen from the fact that I s has only contributions proportional to ϕ e m 4 i with integer m, such that its first Fourier coefficient, and therefore its net drift, vanishes. When a is lowered, in both theory and simulations, attractors form near the domain corners and move inward, since wave length Λ decreases with lower a.
For the cases = = a a 0.06, 0.05, our theoretical prediction does not reproduce the corner attractors, but instead produces a limit cycle close to the domain boundary. Note, however, that our theoretical approximation is expected to break down near the corners of the domain, since we only used an approximation (39) to the exact mechanical Green's function.

Grid pattern for equilibria in a square domain
From figures 6 and 7 we can see that in the square domain, the equilibrium positions the spiral wave's rotation center to occur close to a regular lattice, with spacing equal to half the spiral wavelength, Λ. This situation is similar to the global feedback case described in [48].
From our theory, this regular pattern can be explained as follows. First, due to symmetry, the square's center is stationary. For a given rotation phase, say ϕ = 0, a spiral wave with its rotation center in the middle of the square induces the largest stretch in the region where its tail meets the domain boundary. Suppose this happens at the left side of the square. Now, suppose we move the spiral wave, Λ 2, to the right. Then, the former region of maximal stretch becomes a region of minimal stretch. Since V x is of the form ∫ ϕ ϕ ϕ π f cos ( )d 0 2 , it is mostly affected by rotation phases π 0 and , such that . Conversely, V y is a sine integral and is therefore most sensitive to rotation phases π ± 2, such that . Generalizing to taking m steps of Λ x n 2 in and steps in y, we find m n i( ) Therefore, we expect saddle points on the Cartesian grid of spacing Λ + m n 2 when ( ) is odd. Points with + = m n mod 4 0 will have the same stability as the square's center, and equilibria with + = m n mod 4 2 will have opposite stability. These predictions are confirmed by the dRDM simulation results in figures 6 and 7.
However, in the theoretical prediction in figure 6, the equilibria are not always close to = X Y ( , ) Λ Λ m n ( 2,2). This may be due to approximating the response functions as a delta-distribution, which was performed to produce these plots.  8(b), additional equilibria arise close to the domain boundaries due to saddle-node (SN) bifurcations. Since the SN bifurcation at = a 0.063 interrupts a limit cycle, it can be identified as an infinite-period (IP) bifurcation.

Scaling in large domains
To estimate whether MEF-induced drift is relevant in a practical situation, it is useful to investigate the scaling of the induced spiral wave drift velocity with the domain size, L. Since in the middle of a square domain the velocity vanishes by symmetry, we will try to estimate the drift speed at its local maximum closest to the square's midpoint. From section 6.3, we expect that this point will be found at a position Λ Λ ± ± ( 4, 4)relative to the square's midpoint. The effect of a larger domain size is that the mechanical Green's function needs to be evaluated at a larger distance. Assuming that the dominant contribution of the half-plane mechanical Green's function (38) occurs at = x x 0 and the distance to the spiral core is L, the half-plane Green's function reduces to p L 4 2 , such that the contribution of the four edges is p L 2 . Therefore, the maximal drift speed near the square's center should scale as  We have verified this expression in numerical simulations with various sizes of the square domain ( ∈ L 2 {36, 48, 60, 72}), without changing spatial resolution. Thereafter, drift velocity was interpolated on a grid of size HX = 0.3, and the local maximum closest to the square's center was sought. The outcome is presented in figure 9, which shows that v max indeed almost linearly depends on L 1 2 , as anticipated by expression (41).

Discussion
We have derived a theory of spiral wave drift based on a mechanical Green's function, and we found an expression for the average MEF-induced drift at each point of a two-dimensional domain. Note that in our theory, equation (25) also holds for other geometries.The square domains with fixated boundaries is only one example, since for a different geometry, only the mechanical Green's function needs to be adapted.
Our present approach can be extended to anisotropic excitation of cardiac tissue by local rescaling of the spiral wave; anisotropic mechanical properties and different boundary conditions can be incorporated in the mechanical Green's function.  Due to the presence of a mechanical Green's function, the magnitude of spiral drift is not proportional to the first Fourier amplitude of the excited surface area. In A, however, we show that for tissue of small size ( Λ ≲ L 2 ), the spiral drift pattern can be found from the first Fourier component of the excited area. Note, that in many practically important situations, this assumption can be valid. For example, a study by Nanthakumar et al [49] shows that the typical velocity of wave during ventricular fibrillation in the human heart is 0.67 m s −1 , with an average frequency of 6.8 Hz. This yields a wavelength of 9.9 mm (i.e., about the size of the ventricles of the human heart). In addition, the excited area can be estimated from the experimental recordings using optical mapping experiments [50,51]. Therefore, we hope that our approach, which is based on the first Fourier component, can be applied to experimental data. However, for its application in an anatomically accurate setup, it needs to be extended to three dimensions. Still, we expect from preliminary results that a two-dimensional effective theory taking into account tissue thickness may also be feasible.
Our theory identifies the shape of the wavefront as the main determinant for MEF-induced spiral wave drift. As the shape of the wavefront of a spiral wave is simply given by an Archimedean spiral in an isotropic medium and a rescaled spiral in the presence of anisotropy [52], our theory suggests that these results on spiral wave drift are general. This is also confirmed in our numerical dRDM simulations. It would be interesting, however, to verify the theory using other numerical approaches and models.
Previous analytical work [53] on MEF in cardiac tissue was limited to one-dimensional tissue strands. Here, we developed a theory for two-dimensional media and studied the drift of spiral waves. In two dimensions, the resulting drift field for spiral waves is qualitatively similar to the uniform global feedback chemical system studied in [48]. For example, in a square domain, equilibria of both systems are found on a lattice of spacing Λ 2. Furthermore, works by Zykov et al [54][55][56] revealed that the velocity fields for the spiral wave drift in reaction-diffusion systems in many cases have atypical phase portraits. For example, in [56], unusual equilibrium manifolds of attracting lines have been observed for a system with two-point feedback control. It would be interesting to investigate whether the dynamics reported in these papers can also be found by further variation of the parameters in our system.
Even though our theory is based on linear elasticity, it suitably explains the results of dRDM simulations, which use finite deformations. An extension of our current work could be to find a similar description for nonlinear elasticity, where a crucial step would be to analytically evaluate the deformation field in that case. In addition, an extension of our approach to more realistic settings is straightforward. For that, one can use the approach proposed in [20], which combines an ionic model for human ventricular cells (ten Tusscher-Noble-Noble-Panfilov) and the Niederer-Hunter-Smith model for active tension development.
The results of our study can be tested in an experimental setup. This can be done in micro patterned cardiac cell cultures on elastic membranes [57]. However, such a system may require additional modeling, in which one must consider specific properties of the experimental system, such as the existence of an elastic layer beneath the layer of cardiac myocytes.
An important result of our work is that mechanical feedback and spiral wavelength alone determine the attractors where spirals tend to go. Finding the attractors of spiral waves has become very important as recent cardiac research [58] has shown that ablation of regions where the spiral's tip is located treats atrial fibrillation. Therefore, it is an important task to gain mechanistic insights underlying spiral wave dynamics. We believe our theoretical approach is an important step in the understanding of mechanically caused spiral wave dynamics, and it may lead to new clinical applications.
However, treating the Green's function as a constant is not generally valid, as one can easily see from the counterexample of a half-infinite medium with a zero-displacement boundary. In such a case, the fraction of instantaneous exited area tends to a constant, and omitting the mechanical Green's function would therefore wrongfully forecast the absence of MEF-induced spiral drift. Nevertheless, we can show that this simple approach holds in a square slab of tissue of a size comparable to the spiral wavelength and with fixated boundaries. Although restricted, these conditions apply to mapping experiments in which a thin patch of cardiac tissue is subtended in a rigid frame.
To find such a simplified description, we evaluate (25) for a spiral wave with rotation center X Y ( , )in a square domain, − × − L L L L [ , ] [ , ]. Under no-displacement boundaries, the Navier-Cauchy equation (4) are self-adjoint, leading to the property of mechanical reciprocity that is reflected in the symmetry of the Green's function: = G X Y x y G x y X Y ( , ; , ) ( , ; , ) S 0 0 S 0 0 . By the method of mirror images, the induced stretch field is centered at the points of the boundary closest to X Y ( , )(i.e., the spiral center's orthogonal projections on the boundary). Those four regions lie at distances ∈ d j ( {1, 2, 3, 4}) j , equal to ± ± L Y r L X o from its core; this happens at rotation phases ϕ ϕ π = + j 2 j 0 . In these regions, the overlap integral (35) will be a bell-shaped curve, ϕ ϕ − q d g ( ) ( ) j j 1 1 , with amplitude function q 1 decaying monotonically with the distance, d j . Therefore, the relative spiral drift velocity field can be approximated from (35) as  shape factor, ′ A T , the thickness of the excited zone may be reduced, such that only the length,  ϕ ( ), of the activation front matters: When the square is small enough to accommodate only a single turn of the spiral wave,  attains its local minima at the same rotation phases, ϕ j , as above, whence it is the arc length of a spiral between its core and the point with radial distance d j , and ϕ g ( ) 2 is a bell-shaped curve. Thus, using the excited surface area to estimate the MEF spiral drift, one obtains . If these functions exhibit a similar dependency, the MEF system can be approximated by a uniform-feedback system. Figure A1 illustrates that this is possible in a small square with fixated boundaries.