An analytical velocity field of spiral tips in reaction–diffusion systems

Spiral waves are ubiquitous in diverse physical, chemical, and biological systems. The tip (phase singularity) of a spiral wave is considered to represent its organizing center. Here, we derive an analytical velocity field of spiral tips based on the variables of a general two-variable reaction–diffusion (RD) equation. From this velocity field, we can predict the velocities of spiral tips at time t as long as the values of the variables are given at that time. Numerical simulations with two-variable RD models are in quantitative agreement with the analytical results. Furthermore, we also demonstrate the velocity field of spiral tips in the Luo–Rudy model for cardiac excitation.


Introduction
Spiral waves represent a very prominent example of self-organization in quite different active distributed systems, including the Belousov-Zhabotinsky (BZ) reaction [1], the catalytic oxidation of CO on platinum [2], aggregations of Dictyostelium discoideum amoebae [3], eye retina [4], and cardiac tissue [5]. In cardiology, such self-sustained spiral wave activities play an essential role in cardiac arrhythmia and fibrillation [6,7].
The spatiotemporal behavior of a spiral wave is quite complex, but some of its features, e.g., rigid rotation, meander, drift, etc, can be well described by the motion of the spiral's tip-phase singularity, a kind of topological defect [8,9]. The simplest motion of a spiral wave is rigid rotation, in which the tip follows a circular orbit. Quasi-periodic motions of spiral tips are also possible and this phenomenon was called meander [10,11]. Close to the bifurcation, meander occurs generally in two types called outward meandering and inward meandering, dependent of whether their tip trajectories form a flower-like orbit with petals pointing outward or inward, respectively. It is shown that the meandering dynamics of spiral waves is organized around the codimension-2 bifurcation where Hopf eigenmodes interact with eigenmodes resulting from Euclidean symmetry [12]. Experiments on the BZ reaction unfold the bifurcation from simple rotating spirals to meandering spirals in the neighborhood of a codimension-2 point [13]. Under certain conditions, spiral waves can break up due to the strong meandering of their tips [14][15][16].
If the velocity of the spiral tip is known, then its movement can be quantitatively described. Generally, this velocity could be calculated by the locations of the spiral tip at time t and at a subsequent time t + Δt, where Δt T and T is the rotation period of the spiral wave. Many studies have demonstrated that many important features (e.g. meander, breakup, etc) in spiral wave dynamics can be reproduced by two-variable reaction-diffusion (RD) systems [17,18]. This paper uses a general two-variable RD equation and the velocity field of topological defects [19] to derive an analytical velocity field of spiral tips that is independent of the specific models. The analytical velocity field can predict the velocities of spiral tips at time t as long as the values of the variables are given at that time, so it is no longer necessary to have the values of the variables at two moments. Besides simple cases, e.g., the rigidly rotating spiral wave and the meandering spiral wave, the analytical velocity field can also give the velocities of spiral tips in complex cases, such as in spiral turbulence. In addition, the analytical velocity field of spiral tips can be applied to the Luo-Rudy model, a more realistic model of cardiac excitation.

Theoretical results
Theoretical investigations of spiral waves have been performed predominantly in RD systems [20][21][22]. A general two-variable RD system in two dimensions can be described by the following partial differential equations [15,[23][24][25]: where the variables u and v can be chemical concentrations in a hypothetical chemical reactions, or membrane potential and current in a hypothetical physiological medium. The functions f(u, v) and g(u, v) model the local dynamics, e.g. chemical reaction kinetics. D u and D v are the diffusion coefficients of u and v, respectively, The spiral tip can be defined as the intersections of isolines u = u * and v = v * [11,26,27]. Namely, is the locus of the spiral tip. The solution of equation (2) can be generally expressed as Note that, different sets of u * and v * will give different spiral tip locations, and thus can produce different tip trajectories depending on the choice of u * and v * [28]. A phase variable φ can be computed at each point according to [28,29]: The spiral tip can then be described according to the topological charge [30,31], which is defined as where Γ is a closed curve encircling the spiral tip. The medium in each phase of the activation-recovery cycle encircles the spiral tip, but at the tip itself, the phase is undefined, and the tip is thus a phase singularity. A phase singularity appears when the line integral of the phase change around a point is 2π or −2π, i.e., n t = 1 or −1.
Topological current theory [19,32] shows that Here, ∂S = Γ is the boundary of the surface S, D 0 (u/x) is the Jacobian determinant, and the δ-function is defined as . So, the charge density of the spiral tip is It can be proved that (more details can be found in the appendix) where V 1 and V 2 are calculated according to equation (7). The theoretical descriptions oscillate around the numerical results due to the discrete error. Model parameters are a = 0.84, b = 0.07, and ε = 0.04.
where D 0 (u/x) tip is the value of D 0 (u/x) at the spiral tip.
From topological current theory, we can also obtain the continuity equation satisfied by the charge density (more details can be found in the appendix): where ∇ = i ∂/∂x+ j ∂/∂y, and V(x, y, t) is the velocity field of spiral tips: Note that the above formulas are also applicable to the velocity field of other topological defects, e.g., topological defects in phase-ordering systems [33] and wave dislocations in light beams [34]. In equation (6), to calculate ∂ t u and ∂ t v, we must know u and v at time t and at a subsequent time t + Δt (Δt T). Interestingly, from the general two-variable RD system (1), we can calculate ∂ t u and ∂ t v from u and v at time t. So, we can predict the velocities of the spiral tips at time t as long as we know the values of u and v at that time: The expressions given by equation (7) for the velocity field of spiral tips are useful, because they avoid the need to explicitly specify the positions of the tips.

Numerical simulations
To test the validity of the velocity field of spiral tips in equation (7), we use the Bär model [15] (its supporting spirals have wealthy behaviors including rigid rotation, meander, breakup, etc) to numerically study four different types of spiral wave dynamics: rigidly rotating spiral wave, meandering spiral wave, generation of a spiral pair, and annihilation of a spiral pair. The model kinetics is as follows: The simulations are conducted on 1024 × 1024 grid points employing the explicit Euler method. The space and the time steps are Δx = Δy = 0.025 and Δt = (Δx) 2 /5, respectively. No-flux conditions are imposed at the boundaries. In this paper, to ignore the boundary effect, the regions near the boundaries are not considered. Figure 1(a) shows the u field of a rigidly rotating spiral wave at a certain time, and figure 1(b) shows the field of D 0 (u/x). We can see that the Jacobian determinant D 0 (u/x) is localized at the rotation center, as D 0 (u/x) has a maximum value approximately at the rotation center and is almost zero in other regions [35]. The isocontours of the two-state variables intersect only near the rotation center, and are parallel to each other far from the rotation center. That is, when far away from the rotation center, equation (2) has no solutions, then according to the theorem of implicit function, D 0 (u/x) = 0. Thus, tips only exist near the rotation center. Due to the discrete error in the numerical simulations, D 0 (u/x) can be very small but not equal to zero in the areas far from the rotation center. Therefore, we limit the absolute value of D 0 (u/x) to be greater than a certain value (can be determined by the numerical simulations) when calculating the values of V 1 and V 2 . We find that all tips (different sets of u * and v * ) rotate rigidly.
We calculate the velocity field according to equation (7), and figure 1(c) shows the field distributions of velocity V. Since spiral tips only exist near the rotation center, we can see that the velocity vectors V only exist around the rotation center of the spiral and do not occur in other regions, and that the velocity field is rotationally symmetric about the point where V = 0. Figure 1(d) shows three different isolines of u * (solid lines) and three different isolines of v * (dashed lines). The u * isolines intersect the v * isolines and produce nine intersections that can all be regarded as tips of a spiral.
Additionally, from equation (7), we calculate V and determine the location where V = 0 (marked as '•' in figure 1(e)). Figure 1(e) shows the tip orbit (red circle) of the spiral and the orbit center (marked as '+'). The location where the velocity V = 0 and the center of the tip orbit is coincident. This result follows from the rotational symmetry. To further confirm the validity of equation (7), figure 1(f) compares the values of speed V obtained by theory with those from direct numerical simulation.
We consider a meandering spiral wave in figure 2(a). Similar to the case of the rigidly rotating spiral wave, D 0 (u/x) is localized, and there is a rotating velocity field for the case of the meandering spiral wave, as shown in figures 2(b) and (c). Interestingly, there also exists a point where V = 0 in the velocity field, which is moving in space with time. The velocity field is not rotationally symmetric about the point where V = 0, which is different from the case of the rigid rotation (see figure 1(c)). Figure 2(d) shows the meandering tip trajectory (red curve) with outward petals [12], which is given by the intersection of the isolines u * = 0.50 and v * = 0.35. The trajectory of the point where the velocity V = 0 is shown by the blue curve. Compared with the spiral tip, the point V = 0 can also be considered to represent the organizing center of the spiral wave. Now, we consider a spiral wave which is breaking up [18], as shown in figures 3(a) and (b). Figure 3(c) shows that, when a new spiral pair is generated, the velocity directions near the breakup region (i.e., the area enclosed by the dashed rectangular) are opposite and the speeds are very large. The open ends will become the new organizing centers for spirals later in the simulation. These findings agree with the theoretical results of figure 6(a) in [19], which show that the speeds of topological defects are very large when they are generated.
In addition, we can also find some regularities in the breakup process by observing the distributions of D 0 (u/x), as shown in figure 3(d). Specifically, the maximum of D 0 (u/x) is located at the original spiral tip. When new spirals are generated (the upper part in the panels), the value of D 0 (u/x) is positive on one side and negative on the other side, which means that the topological charges of the two new generated spirals are opposite due to n t = sgn D 0 (u/x) tip [35]. Furthermore, the value of D 0 (u/x) in the intermediate region between the two sides equals zero; this finding is consistent with the results from the branch process of the limit point in section 3.3 in reference [19].
To better observe the velocity vector distributions, we magnify the dashed rectangular areas of figures 3(a) and (c), as shown in figure 4. The left, the middle and the right panels show the distributions of the u field, the velocity field and the isolines of u and v, respectively. From the patterns of the u field in  figures 4(a) and (d), it seems that there are no tips, but through the corresponding isolines of u and v as shown in figures 4(c) and (f), we can see the intersections generated by the isolines of u and v. That is, tips exist in these regions.
In figure 5, we also study the case of an annihilation of a spiral pair. Figures 5(a) and (b) show the evolution of the u field for a spiral pair. Figure 5(c) shows the isolines of u * = 0.50 (blue curve) and  v * = 0.35 (red curve), as well as the velocity vectors. It is found that when the two counter-rotating spiral waves collide with each other and become annihilated, the tip speeds around the annihilating region are also very large, which agrees with the results of figure 6(b) in referene [19]. Also, the value of D 0 (u/x) is positive on one side and negative on the other side, and the intermediate region between them is zero, as shown in figure 5(d). These results are also consistent with the theoretical results from the branch process of the limit point [19].

Discussion
In figures 1-5, we use the Bär model with D v = 0 to demonstrate the velocity field of spiral tips. To check whether the velocity field works for D v = 0, we also consider the case of D v = 0.6 to do the simulations. Compared with figure 1 (D v = 0), similar results are also observed with D v = 0.6, as shown in figure 6.
To further verify that the analytical formula (7) for the velocity field of spiral tips can be applied to different models, we test the Barkley model [25] the FitzHugh-Nagumo model [23,24] , as shown in figures 7 and 8, respectively. It can be seen that there is a rotating velocity field around the rotation center of the spiral, and the velocities obtained by numerical simulations are in agreement with the analytical results.
In order to check whether the analytical velocity field can be extended to more realistic RD models (e.g., ionic models for cardiac excitation), we analyze spiral waves generated by the Luo-Rudy model [36]: where V m is the membrane potential; C m = 1μF cm −2 is the membrane capacitance; D = 0.001 cm 2 ms −1 is the diffusion current coefficient; I ion is the total ionic current, and I ion = I Na + I Si + I K + I K1 + I Kp + I b with I Na = G Na m 3 hj (V m − E Na ) being the fast inward Na + current; is the time-independent K + current; I Kp = G Kp K p V m − E Kp is the plateau K + current; and is the total background current. m, h, j, d, f, and x are gating variables satisfying the following differential equation: where y represents the gating variables. The parameters are simplified as in reference [37]. The Luo-Rudy model is integrated on a 6 cm × 6 cm medium with no-flux boundary conditions via the Euler method. The space and the time step are Δx = 0.0075 cm, Δy = 0.0075 cm, and Δt = 0.00125 ms, respectively. A rigidly rotating spiral is considered as shown in figure 9. To obtain the velocity field distributions, we arbitrarily choose two time-dependent variables in the Luo-Rudy model to replace the variables u and v in equation (7). For example, figure 9(b) shows the velocity field of spiral tips by choosing the membrane potential V m and the gating variable m; figure 9(c) displays the velocity field by using the membrane potential V m and the gating variable x; figure 9(d) exhibits the velocity field by considering the gating variable m and the gating variable x. From figures 9(b)-(d), we can see that there exists a rotating velocity field around the rotation center of the spiral, regardless of the selected time-dependent variables. These velocity field distributions in figures 9(b)-(d) are similar to those in figure 1(c).
Recently, in reference [38], the authors studied the creation and annihilation of a spiral pair in a state of spiral turbulence. They also found numerically that the tip speeds are large when two new spirals are generated or two spirals are annihilated, which is consistent with our results in figures 3-5 and with the topological current bifurcation theory [19]. The generation and annihilation of other topological defects have been studied before. Bray [33] studied the velocity distribution of topological defects in phase-ordering systems and provided a scaling argument associated with the final annihilation of defects, which leads to a large velocity tail. The authors in reference [34] showed that the speeds of wave dislocations in light beams are very large when they are annihilated or generated.

Conclusion
In conclusion, we proposed a quantitative description of the velocity field of spiral tips, providing analytical results that contribute to a deeper and more comprehensive understanding of the dynamics of spiral waves. From a general two-variable RD equation and a velocity field of topological defects, we obtained an analytical velocity field of spiral tips, which is independent of specific models. Using this velocity field, we can obtain the velocities of spiral tips at time t as long as we have the values of u and v at that time. The predictions of the velocity field are in agreement with the results obtained from numerical simulations of RD models and with the topological current bifurcation theory.
Namely, the continuity equation satisfied by the charge density is