Excitable media represent a broad class of non-equilibrium reaction-diffusion systems that play an important role in physical, chemical, and biological applications [1–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–17]. The investigation of spiral wave instability attracts a great attention from a theoretical point of view [18–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.

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₀, v₀). 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 five-point 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.

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 Figures 3, 4.

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.

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 self-sustained 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 in a wider parameter space are shown in Figure 7. The regions of subexcitability (SE) and bistability (BI) are indicated here. Self-sustaining 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.

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.

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

VZ performed the computations. Both authors conceived of the presented idea, discussed the results, and contributed to the final manuscript.

This work was supported by the Max Planck Society and the German Center for Cardiovascular Research (DZHK).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.