Two Domains of Meandering Spiral Waves in a Modified Barkley Model

The stability of rigidly rotating spiral waves is a very important topic in the study of nonlinear reaction-diffusion media. Computer experiments carried out with a slightly modified Barkley model showed that, in addition to one region of instability observed earlier in the original Barkley model, there is another one exhibiting completely different properties. The wave instability in the second region is not related to the Hopf bifurcation. Moreover, hysteresis effects are observed at the boundary of the region. This means that in the vicinity of this region of instability, direct integration of the model equations leads either to a rigidly rotating or meandering spiral, depending on the initial conditions.


INTRODUCTION
Excitable media represent a broad class of non-equilibrium reaction-diffusion systems that play an important role in physical, chemical, and biological applications [1][2][3][4]. For example, wave processes in excitable media are intensively studied in various distributed systems, including the colonies of Dictyostelium discoideum [5], the Belousov-Zhabotinsky chemical reaction [6], the heart muscle [7], the eye retina [8], the neocortex [9], CO oxidation on the platinum single crystal surface [10], and many others.
An excitable medium can be viewed as an ensemble of active elements coupled locally by diffusion-like transport processes. Each individual active element has a resting state, resistant to small external perturbations. However, it can be excited by the application of a suprathreshold stimulus or by interacting with their neighbors. Therefore, locally induced excitation can propagate through the medium as a self-sustaining wave. Such a wave represents a rapid transition from a stable resting state to an excited one followed by a slow recovery transition (refractory) back to the resting state. Under normal conditions, the wave back follows the wavefront, and they never touch each other.
However, under some special conditions, the propagating wavefront can be broken [1,11]. Then the front and the back of the wave propagating in a two-dimensional medium coincide at one point called a phase change point [2]. Near this point, the front and the back are moving in opposite directions and the boundary of the excited region curls around this singularity point. As a result, the broken wave is winding up into a spiral permanently rotating within the medium.
Such self-sustained activity unexpectedly appearing in cardiac or neuronal tissues strongly destroys their dynamics that results in life-threating diseases. In this context, an understanding of possible scenarios of spiral wave dynamics is of great theoretical importance and has many practical applications.
One important aspect of this study is investigation of spiral wave stability. In a homogeneous low excitable two-dimensional medium spiral wave rigidly rotates around a round core. However, under a variation of the medium's parameters this well-ordered dynamics can be destroyed that leads to a transformation of a circular trajectory of the spiral wave tip into the so-called meandering one, e.g., hypotrochoid or epitrochoid [12,13].
Spiral wave meander has been observed in experiments with chemical solutions [14] and in computations performed with different reaction-diffusion models [15][16][17]. The investigation of spiral wave instability attracts a great attention from a theoretical point of view [18][19][20][21].
In this study we would like to find out domains of the spiral wave meandering within the parameter space of a slightly modified Barkley model of an excitable medium.

MODEL
In many studies it was demonstrated that the basic features of the wave dynamics can be reproduced by the two-component reaction-diffusion system where the variables u and v represent the activator and inhibiter species, respectively. Typically the nullcline F(u, v) = 0 is a non-monotonic function creating possibility for undamped wave propagation. The second nullcline G(u, v) = 0 is monotone and intersects the first one at only one point (u 0 , v 0 ). Below the functions F(u, v) and G(u, v) are taken in the form proposed by Barkley [22]: Note, that in the original Barkley model the value of the parameter k ǫ is fixed as k ǫ = 1. Three other constants a, b, and ǫ have been used as important control parameters. A variation of each of these three parameters results in a simultaneous influence on such important medium's characteristics as the propagation velocity, pulse duration and refractoriness. In the modified model under consideration the constant k ǫ is introduced, which has no influence on the duration of a single pulse and its propagation velocity. However, this parameter allow us to control the recovery process because its characteristic time is determined as the product k ǫ ǫ. Thus the activation and the recovery processes have different time constants, if k ǫ = 1. Such a jump in the characteristic time constant is a fairly common and useful tool in simulations of excitable media [12,23,24]. In all computations below the parameter D is fixed as D = 1. The Laplacian in Equation (1) was approximated using the fivepoint finite-difference method on the rectangular 500 × 500 grid with spatial step x = y = 0.3. After this spatial discretization the model equations are integrated in time with the explicit forward Euler method with time step t = 0.01 and no-flux boundary conditions. The spiral wave tip is determined as a point where u = 0.5 and du/dt = 0. A part of an isoconcentration line u(x, y, t) = 0.5 corresponds to the wave front where du/dt > 0, and another part, where du/dt < 0, represents the wave back.

Single Domain of Meandering Spiral Waves
As the first step of our study the parameters are fixed as ǫ = 0.01 and k ǫ = 2, while the constants a and b are used as important control parameters. The obtained computational results are illustrated in Figures 1, 2. Within the parameter space shown in Figure 1 there is a line which determines the boundary of the bistability (BI) domain, where b < a − 1. Here the nullclines of the system (1)-(3) have two intersections points. We analyse another part of the parameter space, where b > a − 1 and the system has only one rest point. Here he existence of spiral waves is limited by another line, where the radius of the core of the spiral wave becomes infinitely large. In the subexcitable (SE) region above this line, the wave segments formed after a wave break are not able to curl around the created singularity point, but simply shrink and disappear. An analytical expression for this line was obtained earlier [25] and has the form where , as it was shown in [26]. Here d u and c p are the duration and the propagation velocity of a single pulse through a one dimensional medium, correspondingly. It can be seen, that the analytical approximation expressed by Equation (4) is in good agreement with the direct reaction-diffusion calculations illustrated by asterisk in Figure 1.
In order to analyse the dynamics of the spiral wave, numerous calculations were performed at various points in the parameter space. A broken plane excitation wave [2] was used as initial conditions. Initially, we fixed a relatively small value of the parameter b. A rigidly rotating spiral wave with a large core was generated near the boundary of the SE region. Then the parameter a increases step by step from one calculation to the next. The size of the core decreases as a increases, and the rotation period decreases. At some computational step, rigid rotation becomes impossible, and a meandering trajectory of the spiral wave tip is observed. This occurs on the left boundary of the light gray region in the Figure 1.
Meandering spirals were observed in the entire light gray region. It is found that in this meandering region the trajectory of the spiral wave tip may look like an epitrochoid (Figure 2A) or a hypotrochoid (Figure 2B). In the white region, to the right of the light gray region and until the BI domain, the tip of the spiral wave moves along a circular trajectory. The radius of this trajectory strongly decreases as a increases.
The computational data shown in Figures 1, 2 look qualitatively similar to ones obtained earlier for the original Barkley model with k ǫ = 1 and ǫ = 0.02 [22,27]. However, the size of the instability domain is considerably smaller in the case under consideration. Note, that while the used value of the parameter ǫ is smaller, the characteristic recovery time determined by the product k ǫ ǫ remains the same.

SECOND DOMAIN OF MEANDERING SPIRAL WAVES
In the second part of our study the value of the parameter ǫ is further decreased to ǫ = 0.005 and k ǫ is increased to k ǫ = 4 in order to conserve the characteristic recovery time. The data obtained in the corresponding computations are shown in  As well as in the previous case, the existence of spiral waves here is limited by the lines defined by Equations (4) and (5). Note that the accuracy of the analytical estimate represented by Equation (5) becomes better as ǫ decreases. Figure 3 clearly shows that there are two regions of instability in the parameter space. The spiral waves in the light gray region show dynamics very similar to those observed in the light gray region in Figure 1. Here, the trajectories of the spiral wave tips resemble epitrochoids or hypotrochoids, for example (see Figure 4A).
In the dark gray region, the tip trajectories become much more complex and are not as well ordered as shown in Figure 4D. In the parameter region surrounding these two regions, the trajectory of the spiral tip is circular.
Note that the light gray domain in Figure 3 is much smaller than in Figure 1. You can also see that the radius of the circular trajectory of the spiral tip is much smaller for these ǫ and k ǫ values, while the values of a and b are the same. This follows from a comparison of Figures 1B,C with Figures 3B,C.

HYSTERESIS PHENOMENON
As the next step of our study the value of the parameter ǫ is considerably decreased to ǫ = 0.001 in the numerical computations. Simultaneously the parameter k ǫ is increased to k ǫ = 20 in order to conserve the characteristic recovery time.
Under these modified values a part of the parameter space shown in Figure 5 looks qualitatively similar to the picture obtained for the original Barkley model as well as for one shown in Figure 1. Within the subexcitable region SE there are no selfsustained spiral waves. Wave segments initiated in this parameter region are shrinking and disappear. Within the rest of the parameter space presented in Figure 5 self-sustained spiral waves have been observed. They are rotating rigidly within the white region, while inside the light gray region they are meandering. Some examples of spiral wave dynamics are shown in Figure 6.
However, this is only a very small part of the entire parameter space investigated at these parameter values. The results obtained  Within the white region between these two domains rigidly rotating spirals with a circular core have been observed. Black spots correspond to parameter values used in Figure 8.
in a wider parameter space are shown in Figure 7. The regions of subexcitability (SE) and bistability (BI) are indicated here. Selfsustaining spiral waves are observed between these two regions. Within the narrow white region, the rigid rotation of spiral waves is stable. The transition to meandering spiral motion occurs in a very small light gray region with a ≪ 1 and b ≪ 1, which is almost invisible in Figure 7 but is shown in Figure 5.
In the dark gray region, the trajectories of the spiral tips are very different from those of the hypotrachoids and epitrachoids shown in Figure 6. A step by step increase of the parameter a within the dark gray domain results in a strong transformation of the spiral tip trajectory. Indeed, rigidly rotating spiral describing  a perfect circular shown in Figure 8A transforms into a jagged trajectory in Figure 8B. Further increase of a results in increasing of angular loops of the trajectory in Figure 8C and their dynamics becomes more irregular in Figure 8D.
Moreover, at the boundary of this region a hysteresis effect in the spiral wave dynamics has been observed. This phenomenon is illustrated in Figure 9. Here the trajectories of the spiral wave tip are shown obtained for different values of the parameter a and b = 0.015. The computations have been started at a = 0.31 and result in rigidly rotating spiral shown in Figure 9A. This stationary rotating wave is used as the initial conditions for the next computations performed with a = 0.32 and illustrated in Figure 9B. After a short transient process the spiral wave trajectory approaches the circular shape. However, a jump to a = 0.325 leads to a destabilization of the rigid rotation and appearance of a rather complicated trajectory, shown in Figure 9C. This wave pattern has been used as the initial conditions for the computations in which the parameter a has been returned back to a = 0.32. However, the spiral tip trajectory does not return back to a circular one, as can be seen in Figure 9D. A rigid rotation restores only for a = 0.31. The further decrease of a also results in a circular trajectory. Thus, it is demonstrated that for a = 0.32 the shape of the spiral tip trajectory depends on the initial conditions.
The observed hysteresis effect exists not only for b = 0.015, but for all other values of b corresponding to the boundary of the instability domain represented by a dashed-dotted line in Figure 9. In particular for a = 1.0 and b = 0.328, as well as for a = 1.4 and b = 0.51. It has been observed not only by a variation of the parameter a and fixed parameter b, but also by a variation of the parameter b and fixed a.

SUMMARY
Thus, the numerical computations performed with a slightly modified Barkley model demonstrate the existence of two quite different parameter regions of spiral wave instability. Within a region located near the SE domain a transition from rigid rotation to spiral meandering follows a well known scenario. Here the instability is induced by the Hopf bifurcation that results in a hypotrachoidal or epitrachoidal trajectory of the spiral wave tip.
The spiral tip trajectories look more complex in the new found region (see Figure 4). The smooth circular trajectory transforms here into a jagged one and even becomes randomized (see Figure 8). This resembles a transition to hypermeandering spiral dynamics reported for the FitzHugh-Nagumo model [13], but is very unusual for the well studied Barkley model. The observed instability cannot be explained by the Hopf bifurcation as was done for the original Barkley model.
At the boundary of this new found instability region the hysteresis phenomenon was detected (see Figure 9). Note, that the similar hysteresis phenomenon was recently observed in the context of the Barkley model within the bistability region [25]. Moreover, a hysteresis phenomenon has been described in context of the FitzHugh-Nagumo model [28,29].
Thus, the results obtained are quite general and applicable to quite different reaction-diffusion models, which should stimulate further research in this area.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.