Efficient control of spiral wave location in an excitable medium with localized heterogeneities

We show that a spiral wave core can be guided by feedback control through a two-dimensional (2D) medium along a virtual 1D detector of given shape. To this aim, short perturbations of excitability are applied globally to the medium each time the spiral wave front is tangent to the detector, or touches its open ends. This relatively simple and robust feedback algorithm is realized in experiments with the light-sensitive Belousov–Zhabotinsky (BZ) medium and in numerical simulations of the underlying Oregonator model. A theory is developed that reduces the description of the spiral wave drift to an iterated map from which the drift velocity field for the motion of the spiral core can be obtained. This drift velocity field predicts both the transient as well as the stationary trajectories of the drifting spiral waves in good agreement with experimental and numerical data. It is shown that the drift velocity is limited by instabilities which arise under high perturbation strength or large delay time. We propose a method to suppress the observed instabilities in order to increase the velocity of feedback mediated resonant drift. Our results might be useful for the control of spiral wave location in a wide variety of excitable media.


Introduction
Spiral waves represent famous examples of self-organization in spatially extended nonequilibrium systems. They have been observed in excitable media of quite different nature, including aggregating slime mould cells [1], the Belousov-Zhabotinsky (BZ) reaction [2,3], heart muscle [4,5], retina of eyes [6], CO oxidation on platinum [7] and frog eggs, through which calcium waves are propagating [8]. The controlled motion of spiral waves induced by a modulation of medium excitability is important for many applications, e.g. for the low-voltage defibrillation of cardiac tissue [9]- [11].
Depending on the parameters of a spatially homogeneous excitable medium under stationary conditions, a spiral wave tip describes a meandering or circular trajectory near a motionless core center [12,13]. Spatial or temporal modulation of excitability is actively used in many studies in order to clarify fundamental features of spiral wave dynamics and to create the basis for possible applications [14]- [20]. In particular, resonant drift of spiral waves induced by periodic temporal modulation of excitability has been proven to be a very generic phenomenon in excitable media [21]- [23]. Very promising control strategies rely on feedbackmediated parametric modulation. In this case, the modulation period always coincides exactly with the actual rotation period of the spiral wave. Moreover, while under external periodic forcing the direction of resonant drift depends on the initial orientation of the spiral wave, the direction of feedback-induced resonant drift is solely determined by the location of the spiral wave core. In consequence, using an Archimedean approximation for the spiral wave shape, the underlying reaction-diffusion model supplemented with terms accounting for the feedback loop can be reduced to an autonomous dynamical system for the coordinates of the spiral core. The corresponding velocity field of feedback-mediated resonant drift simplifies the analysis considerably. Attractors of this velocity field depend on experimentally accessible parameters in the control loop such as the feedback strength or the delay time [24,25].
The most elaborated feedback algorithm to control spiral waves is so-called one-channel feedback [9,26,27]. In this case, a short perturbation of excitability is generated each time a spiral wave front passes through a particular measuring point of the medium inducing the drift of the spiral wave core along a circle centered at the measuring point. Here we report experimental and theoretical results for a novel feedback algorithm which assumes that a short perturbation is generated each time a spiral wave front is tangent to a given curve or touches its open end. These results are supplemented and confirmed by numerical simulations of the underlying Oregonator model. In summary, we state the following general problem in this paper: suppose there is a 3 freely rotating, unpinned spiral wave at a certain position in an excitable medium. How can we move the spiral core as fast as possible along a chosen trajectory to any desired position avoiding those areas where it could be captured by defects or heterogeneities?
The active layer is illuminated by a video projector (Liesegang dv560flex) with intensity controlled by a computer. The light is filtered with a band pass filter (BG6, 310-530 nm). Every one second, images of oxidation waves appearing in the gel layer are detected in transmitted light by a CCD camera (Sony AVC D7CE) and digitized with a frame grabber (Data Translation, DT 3155) for immediate processing by the computer. During the same time step the light intensity generated by the projector can be changed in accordance with the feedback algorithm.
A single spiral wave, which constitutes the initial condition for all experiments, is created near the center of the gel disk by breaking a wave front with an intense light spot. For the chosen concentrations and at a steady level of light intensity (0.36 mW cm 2 ) the spiral wavelength is λ 0 ≈ 2.0 mm. Far from the core center the rotation period is T ∞ ≈ 40 s. The location of the spiral wave tip is defined online as the intersection point of contour lines (0.6 × amplitude) extracted from two digitized images taken with time interval 2.0 s. In the absence of feedback the spiral tip describes a four-petal hypocycloid-like trajectory.
The dashed red line in figure 1 represents an example of a virtual detector used in the feedback algorithm under investigation. A spiral wave (not shown in the figure) is initially located to the right of the detector. At each instant a wave front is tangent to the detector a short light pulse (duration 2 s) is generated. Application of this pulse changes the excitability of the medium, modifies the unperturbed trajectory of the spiral tip, and induces a shift of the spiral core. Due to the rotation of the spiral a sequence of such pulses will be generated and all particular shifts will be accumulated in time resulting in a spiral wave drift.
In our experiments a straight line, a line segment and curved segments have been used as detectors. Figure 1 shows three trajectories of the spiral tip observed in the case of a line detector. After a short transient (not shown here) the spiral core is drifting parallel to the detector line towards the boundary of the reactor. Small variations in the initial location of the spiral core result in a slightly different transient process, but finally the drift of the core center occurs along the same path. However, if the initial location of the spiral core is shifted by a distance of about one wavelength, the spiral core drifts along a quite different path. The attractor structure observed in our experiments is clearly seen in figure 1. Although the spiral wave tip describes a rather complicated meandering trajectory consisting of many loops, the average drift of the core center occurs along a straight line parallel to the detector. The distance between neighboring paths is roughly the wavelength λ 0 of the spiral. Thus, the induced drift allows us to push the spiral wave to the boundary of the medium along the shortest pathway.
Moreover the proposed feedback algorithm can be applied also to shift the spiral wave core from a given point A to a desired position at another point B along a curved detector. This might be important in media with localized defects. Indeed, let us assume that the direct connection between A and B is blocked by a defect D, as shown in figure 2. A spiral core moving along a straight line will be anchored at the defect (see figure 2(b)). However, applying a curved detector as illustrated in figure 2(a) we can shift the spiral wave core from A to B, and avoid it becoming captured by the defect. In these experiments short excitability perturbations have been generated with a time delay τ after the instant when the spiral wave front touches the detector. The time delay induced in the feedback loop allows us to reduce the distance between the attractor line and the detector and makes the control much more precise.
In the following section a theoretical description of the observed phenomena is elaborated to analyze the limitations of the proposed control method and to optimize its efficiency.

Drift velocity field near a straight line detector
A very important property of a spiral wave is that its rotation frequency and shape do not depend on how this spiral has been created. The shape of a rigidly rotating spiral can be rather well approximated by an Archimedean spiral as was suggested by Winfree [2]. Then, in a polar coordinate system (r, ) the spiral shape is given as where 0 , ω 0 and λ 0 determine the initial orientation, the angular velocity and the wavelength of the spiral, respectively. Recent computations performed with the Oregonator model [27] and experiments with the BZ reaction [29] also confirm that an Archimedean spiral provides a suitable approximation of the wave front except for a relatively small region of radius r A λ 0 near the rotation center. Moreover, even the shape of a slightly meandering spiral wave exhibits only small oscillations near an Archimedean form and the amplitude of these oscillations vanishes very quickly with r [30]. Therefore, the Archimedean spiral approximation will be used below to obtain the drift velocity field which determines the feedback-mediated dynamics of a spiral wave near a straight line detector. Figure 3 illustrates an Archimedean spiral rotating around the point (x i , y i ). Each time a spiral wave front (thick black line) touches the detector (dashed red line) tangentially, a short control pulse of amplitude A is applied globally to the whole medium immediately or with some time delay τ .
The perturbation leads to a shift of the core center. Its new location is given by where h is the magnitude of the displacement and the angle γ i specifies its direction. Due to rotational symmetry, the direction of the induced shift depends on the spiral orientation at the instant the pulse is applied. The spiral orientation is determined by the angle α i plus an additional angle ω 0 τ , which results from spiral rotation with a velocity ω 0 during the time delay τ . The angle ψ between the orientation of the spiral wave and the direction of the induced displacement is a characteristic parameter which depends both on the excitability of the medium and the applied perturbation [27]. Therefore the angle γ i in equations (2) and (3) reads A pure geometrical consideration shows that the angle α i can be expressed as where and The integer k depends on the total delay time τ (defined in equation (23)) in the feedback loop. One ingredient of the total delay is the time delay τ in the feedback loop and another one depends on the spiral location and can be calculated as α(x i )/ω 0 [29].
If we assume that the perturbation pulses produce negligibly small changes of the angular velocity and the spiral shape, the system of equations (2)-(7) represents an iterated map describing the dynamics of the core center. The drift velocity field corresponding to the iterated map (2)-(7) is shown in figure 4.
It is important to mention that according to equations (4)-(7) the direction of the displacement γ i depends only on the distance between the line detector and the spiral core center x i and does not dependent on the coordinate y i . This is consistent with the attractor structure observed experimentally as shown above in figure 1. If the spiral tip is located relatively far away from the detector (x i λ 0 ), than due to r (x i ) ≈ x i the iterated map (2)-(7) can be simplified to the following form The first equation of this map is independent of the second one and exhibits steady states x s satisfying the simple condition This condition can be transformed into an explicit analytic expression for x s where m is a positive integer and T 0 = 2π/ω 0 is the rotation period of the spiral. This steady state corresponds to a resonant attractor which is a line parallel to the straight line detector. In accordance with equation (9), each pulse shifts the spiral center by a distance y = h along this line.

8
In the general case of an arbitrary value of x i one can obtain from equations (2)-(5) the implicit expression for the steady state values x s where r (x s ) is determined by equation (6). Note that this expression can be applied when the Archimedean approximation is valid, i.e. x i > r A . Within this range of x i the difference between x s obtained numerically from equation (12) and from the approximate equation (11) is rather small.
The linear stability analysis of the steady states specified by equation (11) leads to the iterated map for the deviations In general, the value H depends on the spiral location. However, under the approximation Substitution of the ansatz δx i = λ i into equation (13) yields the characteristic equation for the eigenvalues A steady state of the iterated map (8) is stable under the condition |λ| < 1. Generally speaking, the characteristic equation (15) has k + 1 roots, which can be found numerically for a given H . For H = 0 there are k roots λ = 0 and one root is located at λ = 1. For H < 0 this root is located outside the unit circle, indicating an instability. With growing H this root moves inside the unit circle. First, the other roots remain in the vicinity of the point λ = 0. However, under critical conditions two complex conjugated roots cross the unit circle, resulting in the Neimark bifurcation.
The critical value H cr of the dimensionless shift H at the Neimark bifurcation can be obtained analytically. To this end one should take into account that at the bifurcation point the eigenvalue λ can be written as From equation (15) we obtain two equations for the real and the imaginary part Positive values of H cr correspond to the case cos 2k+1 2 β = 0. Substitution of this condition into equation (18) yields The integer l specifies the number of the Neimark bifurcation starting with l = 0 for the first one. The obtained expression for H cr allow us to estimate critical values for the magnitude h of the displacement induced by a single pulse. Approximately, equation (14) gives .
The exact implicit expression of h cr follows from equations (12) and (4)- (7). It reads For any odd values of m the coefficient H determined by equation (14) is negative. This corresponds to unstable steady states separating basins of attraction for stable steady states specified by even values of m (cf figure 4). For 0 < h < h cr (k) these steady states describe the spiral drift along straight lines parallel to the detector as shown in figure 5 (cf figure 1). However, the stability of this linear drift depends on k, which characterizes the total time delay τ . The value of k remains constant within the basin of attraction corresponding to an even value of m. In accordance with equation (11) the explicit expression for k reads where [Z ] denotes the largest integer smaller than Z . For an arbitrary chosen h > 0 the stability of a steady state is broken, if k becomes larger than some critical value. This situation is illustrated in figure 5. When h = 3 is fixed, the instability in the linear drift appears first for k = 4. For h = 6 the instability appears earlier, already at k = 2, and the trajectories of the spiral core become more and more complicated for increasing values of k. The predicted instability limits the averaged drift velocity of the spiral wave as illustrated in figure 6. The linear increase of the drift velocity stops at the bifurcation point and is replaced

Double-pulse feedback
To increase the efficiency of the feedback control we applied an additional negative pulse delayed with respect to the positive one. In figure 7 the time dependence of the feedback parameter I (t) is illustrated. First, a single positive pulse is generated as discussed above. In addition, a negative pulse is generated between two positive ones. For example, for p = 0.5, a negative pulse is added after a half of the rotation period of a spiral wave. This negative pulse causes a displacement of the core center practically in the opposite direction in comparison to the displacement induced by the positive pulse. But after half of the rotation period, the orientation of the spiral wave is turned by 180 • . Thus, the negative pulse applied at this instant leads to a displacement in approximately the same direction as that induced by the previous positive pulse.
The iterated map for double-pulse feedback is similar to equations (2) and (3), but it includes additional terms describing the effect of the negative pulse Here β i is the angle of the spiral orientation at the instant when the negative pulse is applied. It is assumed that each negative pulse induces a displacement −h. For a constant angular velocity of the spiral, the angle β i is given by Generally speaking, the value α i is determined by equation (4). To simplify the following analysis we consider the case x i λ 0 that gives r (x i ) = x i and corresponds to the iterated map equations (8), (9). Then equation (24) can be transformed into The iterated map for the y-coordinate of the spiral core center is derived similarly and reads As well as for single pulse feedback the dynamics of the x coordinate does not depend on the dynamics of the y-coordinate. Following equation (27) the analytic expression for the steady states x s reads where m is a positive integer. This expression coincides with equation (11) for p = 0.5. The steady states with even integer m correspond to resonant attractors which are straight lines parallel to the detector. The steady states with odd integer m describe locations of the separatrixes limiting the basins of attraction. During a stationary drift each pulse shifts the spiral center along the attractor line by a distance y which can be determined from equation (28) as Thus, the maximum of the drift velocity is reached for p = 0.5.
As was mentioned above, the drift velocity achieves its maximum at p = 0.5. In this case the characteristic equation for the eigenvalues reads This equation has k + 2 roots, which can be found numerically for a given value of h. For k = 0 the roots can be found analytically. This case corresponds to the stability analysis of the attractor line closest to the detector. For small positive h this attractor is stable, but loses stability at Comparison of this expression with equation (21) shows that the drift instability appears earlier than for single pulse feedback. In figure 8 the averaged drift velocity of the spiral core center in the y-direction versus the absolute value of the displacement h is depicted. The blue curves correspond to double pulse feedback, computed from equations (24) and (25). For comparison the corresponding dependence for single pulse feedback is shown by red lines. It can be seen that for most of the attractors the drift velocity induced by double pulses is larger than by single pulses. The only exception is the case k = 0. Here, under double pulse feedback the instability occurs at a three times lower control strength than under single pulse feedback, and the maximum drift velocity is smaller. However, even in this case the described double pulse algorithm increases the efficiency of the feedback control by application of small modulation pulses.

Feedback-mediated resonant drift observed in reaction-diffusion systems and in experiments
Now, our aim is to compare the results obtained from the iterated map with numerical simulations performed with reaction-diffusion models and with the experimental data.
The light-sensitive BZ reaction applied in our experiments is modeled by the modified Oregonator equations [31] Here the variables u and v correspond to the concentrations of the autocatalytic species HBrO 2 and the oxidized form of the catalyst, respectively. The inhibitor v does not diffuse, because in our experiments the catalyst is fixed within the gel layer. The term I = I (t) corresponds to the the light-induced bromide flow which is assumed to be proportional to the light intensity. The term I (t) is used as control signal in our numerical simulations. All computations are performed by the explicit Euler method on a 400 × 300 array with a grid spacing x = 0.2 and time steps t = 0.002. The autonomous system with I (t) = 0.01 has a steady state which is stable with respect to a small perturbation. However, a supra-threshold perturbation, once locally applied, gives rise to a concentric wave propagating through the medium. A spiral wave rotating counterclockwise near the center of the simulated domain was created by a special choice of initial conditions. The spiral tip performed a compound rotation (meandering motion) including at least two different frequencies. The oscillation period measured far enough from the symmetry center of the unperturbed tip trajectory was T ∞ = 7.65.
In accordance with the first feedback algorithm a short pulse increasing the value I (t) is generated each time the rotating spiral wave front is tangent to a straight line used as a detector. The duration of this pulse was fixed to 0.7 and is small in comparison to the rotation period.
First of all, we performed computations with relatively small pulse amplitude (A = 0.002) in order to check the applicability of the iterated map (2)-(7) to the Oregonator model. Several trajectories corresponding to different initial locations of the spiral center are shown in figure 4. It can be seen that the iterated map computed with ψ = 2.39 predicts quantitatively the asymptotic values of the attractor distances from the detector line. Moreover, the transient parts of the drift pathways are in good agreement with the map predictions. Figure 9(a) shows the spiral tip trajectory obtained numerically under this feedback control for the attractor characterized by k = 1. The stationary drift depicted in this figure is established after a short transient which is not shown here. Thin green curves indicate the trajectories of the spiral wave tip. The time average of this trajectories corresponds to the paths of the core center and is shown by thick black curves. Until the amplitude A of the applied pulses is rather small (A = 0.003), the drift of the core center occurs along a straight line. For larger amplitude (A = 0.0096) the pathway of the core center oscillates near a straight line with relatively small amplitude a (cf figure 5). These oscillations appear due to the instability predicted by the iterated map. Indeed, the magnitude h of the single shift of the core center is a monotonously growing function of the amplitude A. Therefore the instability should appear when A > A cr . The instability becomes more pronounced with growing amplitude of pulses. The pathway of the core center obtained for A = 0.011 deviates very strongly from a straight line. Therefore, the average drift velocity in the y-direction decreases. The average drift velocity induced by the applied feedback control is shown in figure 6 as a function of the magnitude h of a displacement induced by a single pulse. These numerical data are in good quantitative agreement with the results obtained for the iterated map (2)-(7).
Another prominent example of a reaction-diffusion system supporting rotating spiral waves is the FitzHugh-Nagumo medium [32,33] given by the equations Here, u(x, y, t) and v(x, y, t) represent the dimensionless concentrations of the activator and the inhibitor, respectively, β, χ, and 1 are given parameters and D u denotes the diffusion constant. This model is also widely used to simulate cardiac dynamics [34]. In this case u represents the transmembrane potential and v specifies the ionic current. The term F(t) specifies the control signal applied to the medium. Below the model parameters are fixed as = 0.14, β = 1.2, χ = 0.5 and D u = 1.0. The computations are performed by the explicit Euler method on a 1000 × 400 array with a grid spacing x = 0.5 and times steps t = 0.02. The control signal F(t) is generated in accordance with the single-pulse feedback algorithm described above. That means, each instant the spiral wave front touches the detector line a short pulse of amplitude A is generated. Each single pulse induces a displacement of the spiral core center with a magnitude h. Due to the permanent rotation of the spiral wave it is influenced by a sequence of the controlling pulses. Single shifts are accumulated resulting in a drift of the spiral wave core. In full agreement with the theoretical consideration performed above the drift occurs along a straight line parallel to the detector, if the amplitude A is smaller than some critical value. Such a drift is computed for A = 0.06 as shown in figure 9(b). The appearance of the instability is illustrated by computations performed with A = 0.07. Here the pathway of the core center deviates considerably from a straight line. The amplitude of these deviations becomes larger when the pulse amplitude is increased up to A = 0.08.
These computations demonstrate once more that the spiral wave location can be effectively controlled by the proposed feedback algorithm irrespective of the model of excitable medium used in simulations.
The last test of the applicability of the proposed feedback algorithm is performed during experiments with the light-sensitive BZ reaction. In addition to qualitative results shown in figure 1, we have measured the average drift velocity of the spiral core center as a function of the amplitude A of the applied light pulses. These data are depicted in figure 10. The two curves shown here represent results averaged over several experimental series. During one experimental series the pulse amplitude is changed to obtain about ten different values of the drift velocity. Because the transients are very long, it takes about one day to perform one experimental series. The averaged data shown in figure 10 allow us to conclude that the drift velocity can be increased by increasing of the pulse amplitude. There is an instability which appears when the amplitude exceeds some critical value. This instability limits the drift velocity. The double-pulse feedback algorithm increases the efficiency of the feedback control. Note that all these conclusions are in perfect agreement with the theoretical predictions made above and with the numerical analysis of the reaction-diffusion models.

Discussion
In general, given an unpinned spiral wave at a certain position in the medium, the proposed feedback scheme allows us to move the spiral wave core in a well-controlled manner along a chosen trajectory to any desired position avoiding those areas where the spiral wave could be captured by defects or heterogeneities. Ordinary differential equations for the velocity field of feedback-induced resonant drift are derived within the Archimedean approximation for the spiral wave. The obtained velocity fields predict the movement of the spiral core as observed experimentally and in the numerical simulations.
The theoretical approach is based on an analysis of single shifts of the spiral core caused by perturbing pulses in the feedback loop. The observed motion of the spiral core along complicated trajectories is due to a resonant drift induced by the perturbations applied at suitable phases of the spiral rotation. The main advantage of the proposed feedback algorithm is that it is not necessary to determine these phases artificially, as a result of complicated calculations. They are obtained quite naturally as instances when the spiral front touches a virtual detector of a given shape. This simple rule guarantees the robustness of the proposed control algorithm.
We have found good agreement between the theoretically calculated velocity field of resonant drift and the trajectories of the spiral center obtained numerically or observed experimentally. The developed theory allows us to explain the appearance of instabilities of the induced drift which limit the drift velocity. The understanding of these limitations creates opportunities to optimize the feedback control scheme in order to achieve the quickest drift along a given trajectory.
It is important to stress that the theoretical consideration is based on a very general description of an excitable medium and does not use specific features of the experimental or model systems. The close agreement between the theoretical predictions and numerical or experimental data proves that the results obtained are of general validity and can be applied to the dynamics of spiral waves in quite different excitable media including cardiac tissue. While spatially uniform cardiac tissue is clearly highly idealized and realization of the proposed control algorithm is not possible at the current state of implantable cardioverter defibrillator technology, we hope that our study will be useful for future antiarrhythmic pacing devices.