A theory for spiral wave drift induced by ac and polarized electric fields in chemical excitable media

Spiral waves are shown to undergo directional drifts in the presence of ac and polarized electric fields when their frequencies are twice of the spiral frequencies. Here, we propose a quantitative description for the spiral wave drift induced by weak electric fields, and provide the explicit equations for the spiral wave drift speed and direction. Numerical simulations are performed to demonstrate the quantitative agreement with analytical results in both weakly and highly excitable media.

SCientifiC REPORTs | 7: 8657 | DOI: 10.1038/s41598-017-09092- 6 response functions has been presented in an arbitrary model with differentiable right-hand sides 27,41 . By using the response function theory, the drift laws for spiral waves on curved anisotropic surfaces were investigated. And an asymptotic theory that predicts the drift of spiral waves on general curved surfaces with the anisotropic diffusion has been developed 42 . The response function theory has also been used to study the change of the filament tension for scroll waves caused by a circularly polarized electric field 43 .
Spiral waves in the BZ reaction undergo a directional drift when the frequency of an ac electric field is twice that of the spiral frequency 44 . The direction of the spiral drift changes continuously between 0 and 2π when the initial rotation phase of the spiral varies from 0 to π. A kinematical model for spirals subjected to a strong ac electric field was proposed on a phenomenological basis 45 . This kinematical model has succeeded in capturing many aspects of the drift of spiral waves induced by strong ac electric fields. However, it has not been derived from the underlying reaction-diffusion equations. Thus, its parameters need to be adjusted, and cannot be obtained from the underlying reaction-diffusion equation. And the drift of spiral waves induced by weak ac electric fields has been studied from a theoretical point of view 25,46 . But an analytical and quantitative explanation for this mechanism is needed.
Recently, a polarized electric field that possesses chirality was theoretically proposed 47 and has been implemented in the BZ experiment 48 . It allows us to study the response of spiral waves to a chiral electric field. The drift behavior of spiral waves under the influence of a polarized electric field was investigated numerically. An analytical derivation which neglects the deformation of the spiral is consistent with the numerical results qualitatively 47 . However, an analytical derivation which takes account of the deformation for this mechanism is still lacking. Analytical results directly from the reaction diffusion equation are of importance, for a deeper understanding of the mechanism of the drift and a wider application of the spiral wave control.
In this paper, we derive the drift velocities of spiral waves due to weak ac and polarized electric fields using the response function theory, and compare the analytical results with the velocities obtained in direct numerical simulations for both weak and high excitabilities. The detailed forms of electric fields are shown in Table 1.

Results
The theory of wave propagation in excitable media in the presence of an electric field → = E E E ( , ) x y can be described by reaction diffusion equations, , ( , )] T , and  D and  M are constant matrices. In this paper, we use the FitzHugh-Nagumo kinetics 49,50  such that the medium is highly excitable (see Fig. 1(a)). In the weakly excitable one, parameters are chosen as ε = 0.22, β = 0.78, and γ = 0.8 (see Fig. 1(b)). The spiral tip is defined as the intersections of two isolines of u = 0 and v = 0. In the absence of electric fields, a rigidly rotating spiral wave solution to Eq. (1) is in the form of where → R is the center of rotation, Φ is the initial rotation phase, ω > 0 is the angular frequency, and ρ → − → r R ( ) and ϑ → − → r R ( ) are polar coordinates centered at → R . The spiral tip whose position vector is ζ → t ( ) tip rotates in a circle centered at → R . Without loss of generality, we choose the angle between ζ → = t ( 0) tip and the positive x axis as the initial rotation phase Φ, which is measured counterclockwise from the positive x axis. The rotation direction is determined by the chirality σ, i.e. σ = +1 for the counterclockwise rotating spiral and σ = −1 for the clockwise one 51, 52 . Derivation of the drift equations. In the following, without loss of generality, we only consider the drift of the clockwise rotating spiral, which can be easily extended to study the drift of the counterclockwise rotating spiral. It is convenient that the drift velocity of the spiral is analyzed in the system of reference of the spiral, i.e. the system of reference corotating with the spiral's initial phase and angular velocity ω around the spiral's center of rotation. In this system of reference, the polar angle is given by 0 and Φ′ = 0. According to the response function theory 26 u acts on a rigidly rotating spiral by dc electric fields causing translational and rotational shifts. In particular, a solvability condition leads to the following equation for the drift velocity of the spiral: is the complex coordinate of the instant spiral center, and  R t ( ) is the drift velocity. The inner product 〈w, v〉 stands for the scalar product in the functional space in the system of reference of the spiral (Φ′ = 0) and W (1) is one of the response functions of the spiral. Mathematically, it is the eigenfunction of the adjoint linearized operator corresponding to the critical eigenvalue −iω.
To study the directional drift of spirals, it is convenient to average the motion of the spiral over the rotation period. After the central moving average over the spiral wave rotation period, the drift velocity can be expressed as Thus a drift velocity formula in terms of the electric fields, the response functions, the Goldstone modes, and the initial rotation phase of the spiral is given.
For simplicity, we write the resultant as i (1) ( 1) 2 2 such that  R and Θ can be explicitly written. In Table 2, we compute the drift coefficients μ 1 , v 1 , μ 2 , and v 2 for both weak and high excitabilities using the open source software "DXSpiral" 41 .

DC electric fields.
From the general drift velocity formula in Eq. (5), we first discuss the case of dc electric fields. When a dc electric field is applied, only the first component in Eq. (5) contributes to the directional drift of the spiral. And the drift velocity reads as Substitution of Eq. (6) into Eq. (8) gives the drift speed and direction as  for the two excitabilities, respectively. Quantitative differences between theoretical predictions and numerical simulations are relatively small. The relative differences are less than 0.2% for both weak and high excitabilities. Equation (9) shows that both  R and Θ are independent of the initial rotation phase Φ. This is also observed in numerical simulations. Note that the drift equations of spirals by dc electric fields in Eqs (8) and (9) are identical to the results obtained in refs. 28, 37, 38. AC electric fields. The spiral wave undergoes a directional drift in the presence of an ac electric field when ω e = 2ω (see Fig. 2(a) and (b)). In this resonant case, the drift velocity in Eq. (5) becomes Using Eq. (7), one obtains From above equations, we can draw following conclusions about the resonant drift induced by ac electric fields, which are independent of the special models. Firstly, the drift speed  R is independent of the initial rotation phase Φ of the spiral and the initial phase φ e of the ac electric field, which are confirmed by direct numerical simulations in Figs 2(c),(d) and 3(a) and (b). The analytically obtained values  R agree well with the ones obtained by simulations in both weakly and highly excitable media. Secondly, for the fixed φ e , the drift direction Θ changes continuously when the initial rotation phase Φ varies. And Θ is linear in Φ. The change of Θ is twice as much as that of Φ, i.e., ΔΘ = 2ΔΦ. Thus ΔΘ = 2π if ΔΦ = π. This means that the drift direction Θ keeps invariant when the change of Φ is π. The analytical values Θ are quantitatively consistent with the numerical simulations, as shown in Fig. 2(e) and (f). Thirdly, for the fixed Φ, the spiral drifts in different directions when we change φ e , but the change of Θ is equal to that of φ e , i.e., ΔΘ = Δφ e . Numerical simulations are performed to demonstrate the quantitative agreement with the analytically obtained values Θ in both weakly and highly excitable media, as shown in Fig. 3(c) and (d). Note that the numerical and theoretical results of the spiral drift induced by strong ac electric fields in refs. 44 and 45 showed that  R does depend on Φ, and Θ is not linear in Φ. The reason is clear. Since in refs. 44 and 45 E 0 is equal to 2.0, but in our case E 0 is 0.005. According to the response function theory, the analytical results obtained in this paper are valid for weak electric fields. On the other hand, the relation between Θ and φ e was not studied in refs. 44 and 45. Thus, further laboratory works on the spiral drift induced by weak ac electric fields are expected. We note that spiral wave fronts may break above some electric field threshold 53 .
Polarized electric fields. The general drift velocity formula in Eq. (5) can also be applied to the case of polarized electric fields. In the resonant case (ω e = 2ω), the spiral drifts along a straight line, and the drift velocity reads as A significant feature predicted in Eq. (13) is that when the electric field is circularly polarized and its rotation follows that of the spiral (φ xy = 0.5π), the drift speed  R reaches its maximal value. On the contrary, opposite rotation between the spiral and electric field (φ xy = 1.5π) locks the spiral. This prediction is confirmed by numerical simulations, and the analytically obtained values  R are quantitatively consistent with the ones obtained by simulations, as shown in Fig. 4(a-d). And the drift direction Θ changes continuously when Φ or φ e varies. Their relations are ΔΘ = 2ΔΦ and ΔΘ = Δφ e , which are the same results as in ac electric fields. Moreover, the spirals drift in different directions when we change the phase difference φ xy , and the change of the drift direction is half as much as that of the phase difference, i.e., ΔΘ = 0.5Δφ xy . This unexpected prediction is also confirmed by direct numerical simulations, and the analytical values Θ agree well with the numerical ones, as shown in Fig. 4(e,f).
In conclusion, we have studied the drift of spiral waves induced by weak electric fields analytically, since analytical results are of crucial importance for a deeper and more comprehensive understanding to the mechanism of the spiral wave drift. Using the response function theory, we propose a theory of the spiral wave drift due to weak ac and polarized electric fields. Explicit equations for the spiral wave drift speed and direction in terms of Φ, E 0 , φ e , and φ xy are obtained directly from the reaction diffusion equations, which are independent of the special Figure 4. First row: Drifting behaviors of spirals under the influence of a clockwise (φ xy = 0.5π) and a counterclockwise (φ xy = 1.5π) circularly polarized electric fields with ω e = 2ω, φ e = 0, and Φ = 0. Second row: Dependence of theoretical (lines) and numerical (circles) drift speed on φ xy with φ e = 0 and Φ = 0. Third row: Dependence of theoretical (lines) and numerical (circles) drift angle on φ xy with φ e = 0 and Φ = 0. When the drift speed is 0 (φ xy = 1.5π), the drift angle cannot be defined. Left column: Highly excitable medium. Right column: Weakly excitable medium. models and should be of general significance. These analytical results are quantitatively confirmed by numerical simulations in both weakly and highly excitable media. Although ac electric fields 44 and polarized electric fields 48 have been realized in the BZ system, laboratory works on the spiral drift induced by weak ac and polarized electric fields have not been investigated. We hope that the theoretical results in Eqs (11) and (13) can be verified in experiments.

Methods
In direct numerical simulations of the FitzHugh-Nagumo model in Eq. (1), we use an explicit Euler method and no-flux boundary conditions with the space step at Δx = Δy = 0.05 and time step at Δt = 0.0005 for the grids 500 × 500 in Cartesian coordinate system.
We compute V (1) (Φ′ = 0), V (−1) (Φ′ = 0), and W (1) (Φ′ = 0) using the open source software "DXSpiral" 41 . Since "DXSpiral" is dealing with spiral waves on a polar grid in a disk, to get a similar precision with Δx = 0.05 in direct numerical simulations, we use a disk of radius at 12.5 with 250 radial and 124 circumferential grid cells. Note that in the computation, the initial rotation phase of the spiral in the system of reference of the spiral should be set to be 0, i.e. Φ′ = 0. In our computations, we find that is not. This is related to the fact that the drift velocity of the spiral induced by the dc electric field does not depend on the initial rotation phase of the spiral, while the drift velocity induced by the ac electric field depends on the initial rotation phase.