Spiral waves within a bistability parameter region of an excitable medium

Spiral waves are a well-known and intensively studied dynamic phenomenon in excitable media of various types. Most studies have considered an excitable medium with a single stable resting state. However, spiral waves can be maintained in an excitable medium with bistability. Our calculations, performed using the widely used Barkley model, clearly show that spiral waves in the bistability region exhibit unique properties. For example, a spiral wave can either rotate around a core that is in an unexcited state, or the tip of the spiral wave describes a circular trajectory located inside an excited region. The boundaries of the parameter regions with positive and ‘negative’ cores have been defined numerically and analytically evaluated. It is also shown that the creation of a positive or ‘negative’ core may depend on the initial conditions, which leads to hysteresis of spiral waves. In addition, the influence of gradient flow on the dynamics of the spiral wave, which is related to the tension of the scroll wave filaments in a three-dimensional medium, is studied.


Introduction
Wave processes in active media have attracted a great deal of attention in recent decades. For example, wave processes in excitable media are intensively studied in various physical, chemical, and biological systems, including the colonies of Dictyostelium discoideum [1], the Belousov-Zhabotinsky chemical reaction [2], the heart muscle [3], the eye retina [4], the neocortex [5], CO oxidation on the platinum single crystal surface [6], and many others.
An excitable medium can be viewed as an ensemble of active elements coupled locally by diffusion-like transport processes [7]. The simplest example of an active medium is a bistable distributed system that has two stationary states and can be simulated by a one-component reaction-diffusion model [8]. The propagation of trigger waves of advantageous genes through a distributed population [9,10] or a flat flame propagation through a reactive mixture [11] exhibit these properties. These processes must be considered as transitions from one resting state of the system to another, in which the system will remain for ever. Since a secondary excitation is absolutely impossible, one cannot expect self-sustaining spiral waves in such bistable media.
A more complex example of active dynamics is the so-called excitable medium. Here 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 interaction with its neighbors. Therefore, locally induced excitation can propagate through the medium as a self-sustaining wave. After this wave has passed, the medium returns to its resting state. Thus, repeated excitation is possible, and these active media are capable of maintaining selfsustaining rotating spiral waves [7,12].
To simulate an excitable medium, it is necessary to use at least a two-component model [13][14][15][16]. As a rule, this model has only one resting point, which is stable to subthreshold perturbations. However, in a certain range of model parameters, a second stable rest point may arise. In this range of bistability, the excitable medium exhibits quite different properties than the bistable active system described by the one-component model. For example, the existence of spiral waves has been demonstrated experimentally and in computer simulations [17][18][19][20][21][22].
In this study, we would like to analyse the properties of spiral waves within a bistability parameter range of an excitable medium.

Results
In our computations, we use the well-known and widely applied Barkley model [23]. It is the two-component reaction-diffusion system where the variables u and v represent the activator and inhibitor species, respectively. The functions F (u, v) and G(u, v) are taken in the following forms: Here, the nullcline F(u, v) = 0 is a non-monotonic function that allows undamped wave propagation. The second nullcline G(u, v) = 0 is monotonic. The intersections of these nullclines determine the number and locations of the resting points. These are defined by two control parameters a and b, while the small constant and the diffusion coefficient D are fixed below as = 0.02 and D = 1.0. The nullclines of the Barkley model for a = 1.0 and b = 0.2 are shown in figure 1(a). In this case, there is only one intersection (0, 0), which corresponds to the single resting state of the system. An application of a suprathreshold stimulus is able to induce propagation of an excitation wave. For the Barkley model with 1, the wave boundary is very sharp and the inhibitor value is considered constant here, while the activator value jumps between zero and one. The propagation velocity c p of a planar wave boundary can be determined analytically as [24] c where v b is the inhibitor value at the wave boundary. If v b = v + = 0 near the wave front, its velocity in the normal direction is positive c + = √ 2D(a/2 − b)/a. The propagation velocity of the wave back c − should be negative, because it is following the wave front. However, its absolute value must be equal to c + to preserve the stationary waveform. Therefore, the inhibition value v − at the back of the planar wave can be estimated This dynamics of the two variables is shown schematically in figure 1(a) by dashed lines indicating a fast jump from the resting state toward the right branch of the nullcline F(u, v) = 0, then a slow motion along this branch until v − , then a jump to the left branch and a slow return to the resting state. For this set of parameters, there are many methods to generate a rigidly rotating spiral wave, shown in figure 1(b) [25]. The gray region here corresponds to a part of the wave where u(x, y, t) > 0.5. Such a counterclockwise rotating spiral wave is a typical example of a self-sustaining activity in the Barkley model, which has been studied in detail, for example, in [26]. In this work and many other papers, it was postulated that the spiral waves do not exist within a bistability parameter range, which can be expressed analytically as b < a − 1.
There is a very serious reason for this postulate. Indeed, within a bistability parameter range, a suprathreshold perturbation should induce a transition from resting state to the second excited steady state. Thus, this perturbation will produce a trigger wave of this transition, as in a one-component bistable medium, rather than an excitation pulse. Hence, some methods for spiral wave creation, e.g. listed in [25], do not work within the bistability range, because they require a limited duration of an excitation pulse. However, this postulate contradicts to the fact that a periodic sequence of propagating pulses is possible even in the bistability parameter range [27,28]. Moreover, spiral waves have been experimentally observed in a bistable chemical reaction during the CO oxidation on a Pt(110) surface [18,20] and reproduced in numerical simulations [17,18,[20][21][22].
To understand and overcome this contradiction, we consider the dynamics of the Barkley model with a = 1.25 and b = 0.2, i.e. within the bistability range. The corresponding nullclines are shown in figure 1(c). In contrast to figure 1(a), there are three intersection points. Two of these points, namely (0, 0) and (1, 1), are stable, which corresponds to the bistability region. It means that an application of a suprathreshold stimulus to the resting state (0, 0) does not induce a propagation pulse, but only a trigger wave corresponding to a transition to the excited state (1, 1). However, if one performs computations for this set of parameters with the initial conditions corresponding to the medium's state with the spiral wave shown in figure 1(b), a spiral wave with a changed shape and another rotation frequency is approached, as shown in figure 1(d).
The explanation for this fact is quite simple. Namely, the dashed line in figure 1(c) shows that the system leaves the right branch of the nullcline F(u, v) = 0 earlier than it reaches the second stable point (1, 1). A similar dynamic is shown in figure 1(a), where the jump to the left branch occurs before the system reaches the upper point of the right branch. The reason for this is that the absolute value of the back velocity should be the same as the front velocity. Due to equation (5), this velocity corresponds to the inhibitor value v − = 0.85, which is smaller than the top of the right branch.
Thus, within the bistability parameter range the system dynamics very strongly depends on the initial conditions. This difference, observed in our computations, is very important for spiral wave dynamics.
Note that the core radius in figure 1(d) is much smaller than that in figure 1(b). As the parameter a is further increased, the core radius approaches zero, as shown in figure 2(b). With this set of parameters, there are three intersection points shown in figure 2(a). The unstable point is located at (0.5, 0.5), i.e. exactly in the center of the phase portrait. The entire phase portrait is perfectly symmetric. This symmetry is reflected in the spiral waveform, which resembles the famous Yin-Yang symbol. In fact, the shape of the region where u(x, y, t) > 0.5 is identical to that where u(x, y, t) < 0.5. As a consequence of this symmetry, the size of the spiral wave core is equal to zero.
A further increase of the parameter a leads to another important result, shown in figure 2(d). Here, the trajectory of the tip locates within the region where u(x, y, t) > 0.5, which is usually considered to be excited. Therefore, one can speak about a 'negative' core of the spiral wave. This paradoxical situation can be easily explained if the state u = 1 is considered as a resting one, and a decay to u = 0 corresponds to a 'negative' impulse propagating through the medium [27,28]. It can be seen that the shape of the 'negative' spiral wave in figure 2(d) is very similar to the regular 'positive' wave in figure 1(b). Note, that traditionally the wave front is considered as a part of the wave boundary, where du/dt > 0. However, for a 'negative' spiral wave this part of the wave boundary represents the wave back. This is the reason why the wave front, that is shown as a thick solid curve, is concave near the wave tip. The wave front for a 'negative' spiral wave should correspond to a part of the wave boundary, where du/dt < 0.
In order to generalise these unexpected observations, many computations were performed, which were carried out for different a and b. The results of these computations are summarised in figure 3.
For example, computations performed at the same b = 0.2 with a decreasing a show that at a < 1.0 the core radius becomes larger and eventually approaches infinity. This situation is equivalent to the emergence of a so-called critical finger discovered in [29] and later studied in detail [30][31][32][33][34]. However, the similar scenario is observed when the parameter a > 1.6 and is further increasing. In this case the radius of a 'negative' core approaches infinity.
Two examples of a normal and a negative critical finger are shown in figure 4. In both cases the parameter b is fixed as b = 0.7. In figure 4(a), the parameter a corresponds to the left edge of the existence region, where the radius of a spiral wave core approaches infinity. This wave pattern represents a normal critical finger moving to the right with a constant velocity. In contrast, the pattern in figure 4(b) corresponds to the right edge of the existence region, where the radius of a 'negative' core diverges. This pattern can be viewed as a 'negative' critical finger also moving to the right. The shape of the 'negative' region is similar to the shape of the critical finger in figure 4(a).
The repetition of these computations for different values of b makes it possible to determine the boundaries of the region in which the spiral waves exist, indicated in figure 3 by the lines of triangles and diamonds.  [26]. The brown region corresponds to meandering spirals. The dashed line shows the boundary of the bistability region where b < a − 1. Within the red region, spiral waves rotate around a circular core. Within the blue region, spiral waves have a 'negative' core. The dotted line shows the boundary between the positive and 'negative' core regions expressed by equation (7). The thick solid line corresponds to an orthogonal spiral drift computed with equation (9) containing an advection term. The parameter sets used in figures 1 and 2 are indicated by +. Outside the bistability parameter range, rigidly rotating spiral waves exist within the yellow region while meandering waves observed within the brown region [26]. The dashed line represents the boundary of the bistability region expressed by equation (6). Inside the bistability parameter range, the spiral waves exist only within the red and blue areas bounded by these two boundaries. Within the red region, the trajectory of the spiral tip locates outside the exited region, while within the blue region the core is 'negative'. The dotted line shows the boundary between these two regions with a positive and negative core. An analytical expression for this line can be derived after a corresponding analysis of figure 2(a). Indeed, due to equations (2) and (3), the upper point of the right edge of the u-nullcline is (1, a − b). The distance between this point and the right intersection is a − b − 1. Because of the symmetry of the phase portrait, this gap should be equal to the distance between the left intersection point and the bottom point of the left branch, which is equal to b. Therefore, the expression for the boundary between the red and dark blue areas reads b = (a − 1)/2.
The thick black line within the yellow region corresponds to the boundary between the domains where collapsing and expanding scroll rings can be observed in 3D simulation assuming axial symmetry of the rings [26]. In our computations this boundary is continued within the bistability region. For this aim equation (1) is first modified as where x is the radial distance from the scroll ring axis. In the following 2D computations the position of the core center corresponds to the location of the scroll ring filament. An increase (decrease) of the x coordinate of the core center corresponds to an expanding (collapsing) filament [35]. Note, that to determine the boundary between the collapsing and expanding scroll ring filaments equation (8) can be further simplified as where α is a small positive parameter, which determines the intensity of the advection term. Similar models have been used to simulate chemical excitable media in the presence of an electric field [36]. In our computations the parameter α was fixed as α = 0.002. The presence of the advection term results in the spiral wave drift in the positive or in the negative x directions. Simultaneously there is a drift in y direction. A drift in the positive (negative) x direction corresponds to expanding (collapsing) scroll rings in 3D. The results of computations performed for different values of the parameters a and b show that the boundary between the collapsing and expanding scroll rings penetrates into the bistability domain with increasing b, reaches the symmetry line expressed by equation (7) and then returns back towards small b, as shown in figure 3. Outside the region bounded by this curve the scroll rings expanding in 3D, that corresponds to a negative filament tension and results in the Winfree turbulence of scroll waves [37].
It has been shown that the line determined by equation (7) reflects a symmetry in the reaction terms [28]. It is natural to assume that along this line the core size should tend to zero, as it does in figure 2(b). This assumption works perfectly up to a = 2.4 and b = 0.7. However, with this set of parameters, it was not possible to generate a spiral wave with a zero core radius, because the corresponding wave pattern becomes unstable. Depending on the initial conditions, either a rotation around a circular core can be induced, as shown in figure 5(c), or the spiral wave tip describes a circular trajectory within the excited region, as shown in figure 5(h). This unexpected dependence of the spiral wave dynamics on the initial conditions was observed within a narrow interval 2.398 < a < 2.402 and b = 0.7, as shown in figure 5. For a = 2.398 a rigidly rotating spiral wave with a positive core radius has been observed, as shown in figure 5(a). This wave pattern has been used as initial conditions for computations with the increased parameter a = 2.399 (see figure 5(b)). Further step by step increase of the parameter a until a = 2.401 still results in spiral waves with a positive core radius, as shown in figures 5(c) and (d). However, at a = 2.402 a rigid rotation around a positive core becomes unstable, and a transition to a 'negative' core has been observed, as illustrated in figure 5(e).
After stabilisation of this transition the spiral wave continues the rigid rotation with the 'negative' core radius, as shown in figure 5(j). During a step by step decrease of the parameter a until a = 2.399 this dynamics is conserved, see figures 5(i), (h) and (g). However, at a = 2.398 this rotation becomes unstable and a transition to a positive core has been observed, as illustrated in figure 5(f).
In fact, our computations demonstrate a hysteresis phenomenon within this parameter range. Similar instability and pattern formation as a function of initial conditions have been also observed for larger values of a and b along the line expressed by equation (7).
Note that as a and b increase, the left and right boundaries of this range converge and coincide at a point indicated by an asterisk in figure 3. With this set of parameters, two quite different wave patterns can emerge depending on the initial conditions, as shown in figures 6(a) and (b). In fact, the radius of the spiral wave cores shown in figures 5(c) and (d) diverges here. Above this point at the line determined by equation (7), no stationary rotating spiral waves can be generated.
It is important to stress, that the shape of the region, where spiral waves can be observed, can be approximated analytically. For this aim the important dimensionless parameter can be used where d u is the pulse duration. In the case of a positive critical finger analysed in [29] this parameter reaches the critical value B c = 0.535.
In accordance with equation (3) the pulse duration is expressed as Hence, the parameter B can be expressed as Thus, the curve in the parameter space, where B = B c reads This curve is shown in figure 7 as a thick solid and corresponds to a positive critical finger. In order to analyse the negative critical finger, a symmetrical transformation should be performed as u − = 1 − u and v − = 1 − v [28]. After this transformation the propagation velocity of a negative pulse c n p can be expressed as The duration of the negative pulse reads as According to equations (10), (15) and (16), the curve in the parameter space corresponding to the negative critical finger looks like This curve is shown in figure 7 with a thick dotted line and corresponds to the negative critical finger. Both of these curves correlate quite well with the numerical data obtained above. Note that these two curves approach the symmetry line of equation (7) and coincide at the point where Figure 7. Domaine of spiral wave existence determined numerically (triangles and diamonds) and analytically. Thick solid line corresponds to equation (14). Thick dashed line corresponds to equation (17). Dotted lime corresponds to equation (7).
Thus, the analytical approximations expressed by equations (14) and (17) are in good agreement with the direct reaction-diffusion calculations and give a clear description of the domain boundaries on the parameter plane (a, b), where spiral waves can be induced.

Summary
In summary, the direct computations performed with the Barkley model clearly show that the spiral waves exist in a much wider parameter range than it was previously postulated by Alonso et al [26]. It corresponds very well with the experimental observation of rotating spiral waves in CO oxidation on a Pt(110) surface [18,20] and with a numerical analysis of a modified FitzHugh-Nagumo model [21,22].
Our calculations clearly show that in the region of bistability parameters, spiral waves either rotate around the wave core, which is at rest, or the tip trajectory rotates around the 'negative' core, which is in an excited state. The detailed structure of spiral wave modes is shown in figure 3. The analytical expression equation (7) for the boundary between these two parameter ranges is in quantitative agreement with the results of direct numerical calculations. Along most of this line, the core size is zero. However, in the upper right part, where b > 0.7, this solution becomes unstable. Depending on the initial conditions used, for the same set of parameters here we obtain two different rotation patterns with a normal positive or an unusual 'negative' core. Thus, the duality of the spiral wave core and the hysteresis phenomenon in the dynamics of the spiral wave are observed in this parameter range.
The parameter region in which spiral waves can exist is limited by two lines corresponding to positive and negative spiral waves with an infinitely large core, the so-called critical fingers. Two analytical expressions of equations (14) and (17) have been derived, which agree very well with the results of direct numerical calculations, as shown in figure 7.
It was also important to extrapolate the boundary between the expanding and collapsing spiral ring obtained in [26] to the bistability region, where the spiral waves rotate with either a positive or a negative core.
Thus, the calculations carried out with the help of the Barkley model allowed us to discover many important and unknown properties of rotating spiral waves by variation of two parameters a and b. Another important parameter of the Barkley model is , which was set in the above calculations as = 0.02. It has been shown in many studies that this parameter is very important for the dynamics of spiral waves in monostable [38,39] and bistable media [21,22]. In particular, important results concerning the dynamics of single spiral waves in a bistable medium described by the modified FitzHugh-Nagumo model were obtained in [21]. These results do not contradict the data we obtained for the Barkley model and significantly complement them. We are planning more detailed investigation of the influence of on the spiral wave dynamics in the bistability region.
Since the Barkley model is very general and is widely used to reproduce basic properties of spiral wave dynamics, the data presented here should be applicable to various models and excitable media of different nature. We plan to estimate the parameters of the obtained spiral waves using the previously developed kinematic approach [40].
It should be noted that spiral waves are considered as a probable cause of cardiac arrhythmias. It is important to emphasise that in the region of bistability parameters, spiral waves cannot be initiated according to many known scenarios, such as those listed in [25]. In our calculations, we used either some special initial conditions, e.g. a broken planar wave, or the spiral wave was initiated first in a monostable medium and then the medium parameters were changed to enter the bistability region. In chemical systems, this method can be reproduced quite easily, for example, by changing the illumination of the Belousov-Zhabotinsky solution. In cardiac tissue, environmental parameters are also highly depended on many external factors and can change over time, for example, as a consequence of arrhythmia caused by a spiral wave of excitation. Relevant possible scenarios should be analysed using more detailed cardiac tissue models.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).