Another type of attractor that can be found in various examples of excitable media (with no-flux boundary conditions) is that of a single spiral wave. In its simplest form of a rigidly rotating wave, the trajectory of the system corresponds to a periodic orbit. Therefore, at a first glance and from the dynamical systems point of view, the attractor seems to be similar to the one considered in Sect. 3.1.1. However, it is known that spiral waves are characterized by phase singularities (see Sect. 2.​1.​7), a topological property which represents a localized wave source that does not require specific anatomical features, as opposed to the plane wave that requires the periodic medium considered in Sect. 3.1.1. This qualitative difference between the two attractors becomes manifest in the spectrum of Lyapunov exponents by the appearance of two additional Lyapunov exponents that are equal to zero: In Fig. 3.4, the subspace belonging to for a spiral wave in the Barkley model is illustrated with three meaningful directions in tangent space. Due to the degeneracy, Lyapunov stability analysis yields three arbitrary Lyapunov vectors in this subspace. In Fig. 3.4, three orthonormalized linear combinations of these arbitrary Lyapunov vectors are depicted, indicating that the two additional zero exponents are a result of the translational symmetry in the system: If the boundary is far away from the rotational center of the spiral—i.e. the two are electrotonically decoupled—a shift of the wave in either of the two available directions will not decay but represents a neutrally stable perturbation that is picked up by the linear stability analysis. The first Lyapunov vector plotted in Fig. 3.4 (second column) corresponds to a perturbation along the trajectory of the system, i.e. parallel to . This is the neutrally stable direction in tangent space expected for a periodic orbit in every autonomous dynamical system, in this case accelerating the dynamics towards a future state. Consequently, the other two linearly-independent directions indeed correspond to a shift of the spiral wave in two different directions. For example, the perturbation in the third column of Fig. 3.4 accelerates wave parts moving (roughly) to upper left of the medium and decelerates those parts moving to the lower right, effectively shifting the whole spiral towards the upper left.

The periodic orbit of a steadily rotating spiral therefore seems to be embedded in an invariant manifold of total dimension three, where the two directions transverse to the flow represent the independence of the dynamics from position of the spiral tip. Although the finding that spiral waves are subject to a Euclidean symmetry (ignoring boundary effects) seems trivial, note that Lyapunov spectra are computed using a spatially and temporally discretized system. In contrast, the translational, Euclidean invariance of the dynamics is a continuous symmetry. A priori, it is at least questionable whether this continuous symmetry is also present in the discretized system. The analysis shown here confirms that this is indeed the case.

The correctness of the numerical approximation in this particular aspect (with respect to the original PDE) depends on the chosen numerical resolution, as shown in Fig. 3.5. For high resolutions, the upper Lyapunov spectrum becomes independent of the numerical resolution and the structure of the neutrally stable subspace is preserved. For low resolutions, the Euclidean symmetry is broken and the neutrally stable subspace becomes one-dimensional: certain positions of the spiral wave are favored and stabilized by the numerical grid, hence yielding negative Lyapunov exponents for the translational perturbations. In other words, the spiral wave tends to pin to the numerical lattice. This effect has already been noticed by comparing spiral tip trajectories for different numerical resolutions [7]. For stable spiral waves without a visible core (i.e. with a stationary tip), however, this strategy is inappropriate, since a retrospective analysis cannot determine whether or not the tip could potentially have been located at a slightly different position. The plot in Fig. 3.5 also indicates that a low numerical resolution does not have identical effects on all directions in tangent space: in contrast to and , higher-index Lyapunov exponents increase towards lower resolutions, meaning that the corresponding perturbations become less stable. This effect is, however, not universal: Fig. 3.6 indicates that for specific parameters in the Fenton-Karma model, a low numerical resolution has a stabilizing effect on , 7, 8] and prevent [9, 10] spiral wave break-up, which would be associated with one or more Lyapunov exponents being elevated above or reduced below zero, respectively. In particular for investigating control strategies for spiral wave activity in excitable media, Lyapunov stability analysis therefore presents a valuable tool to directly assess the accuracy of the numerical approximation. In this way, the faithful reproduction of stability properties in the numerical realization of the PDE can be ensured without relying on the traditional method of checking the sensitivity of indirect indicators to changes in resolution, such as spiral tip trajectories or periods, conduction velocity, action potential duration and the rate of change of the membrane potential during the upstroke phase [11]. The big advantage of the analysis lies in its implicit check for possible perturbations, whereas all other indicators mentioned above only operate on the attractor itself.

In Sect. 3.1.2, it was shown that Lyapunov stability analysis can detect symmetries in an activation pattern. At a second glance, this is quite surprising for the translational symmetries, since the domain under consideration (a bounded region with no-flux boundary conditions) does not have this symmetry. Consequently, it can be concluded that the part of the activation pattern responsible for the observed zero exponents (and possibly also for other exponents in the upper Lyapunov spectrum) is actually much smaller than the domain size and thus protected from the influence of the boundary.² This does not contradict the extension of the corresponding Lyapunov vectors throughout the whole domain (see Fig. 3.4), as translations and rotations—the latter being equal to a perturbation along the trajectory—nevertheless act globally, as long as the spiral wave can reach the whole domain.

One way to assess the extension of the essential part of the spiral in terms of symmetry and stability would be to investigate the dependence of the spectrum on system size. Here, a different approach is chosen: different numbers of spiral waves are initiated in a medium of constant size, using the Barkley model with parameter set B2 (see Table A.1). Figure 3.7a shows example snapshots of the -variable of the Barkley model for one to twelve spiral waves. The images correspond the state of the medium after an initial transient phase of spiral movement and possible disappearance of phase singularities from the medium. From the plot in Fig. 3.7b, it can be concluded that the specific signature of three zero Lyapunov exponents for a spiral wave must be carried by a relatively small region around the phase singularity. For one to twelve spiral waves, the number of zero Lyapunov exponents corresponds exactly to three times the number of spirals in the domain. This means that all the symmetries are preserved for each spiral individually. More specifically, this rules out the possibility that any of the spiral waves are phase-synchronized or that the position of one phase singularity is constrained by interaction with the surrounding ones (in which case some zero Lyapunov exponents should become negative). As indicated in Fig. 3.7c, this additive behavior of the Lyapunov spectrum does in fact extend to the next (more negative) two pairs of Lyapunov exponents. For the number of spirals and parameters considered here, the entire leading Lyapunov spectrum for spirals is therefore identical to that one expected for copies of a system with one spiral wave. This is remarkable, since different actually represent coexisting attractors in the same system. On the one hand, the behavior of the Lyapunov spectrum reveals that these attractors have a very specific structure. On the other hand, it confirms the intuitive interpretation of spiral waves in excitable media as objects.

To emphasize that the directions in the more and more degenerated subspaces in tangent space can indeed be assigned to individual spiral waves, Fig. 3.8 shows four specific Lyapunov vectors from the subspace with growth rate zero in a three-spiral simulation. The first one corresponds to the perturbation along the trajectory, while the other three indicate that the time- or phase-shifting symmetry is indeed preserved for each spiral wave individually, as if they were (almost) independent systems. As an interesting by-product, these specific vectors yield a decomposition of the medium into domains controlled by the individual spirals. At least for this periodic attractor, Lyapunov stability analysis therefore provides a way to objectively define “domains of influence” of individual spirals in contrast to visual inspection of the dynamics [12].

It should be noted that the complete preservation of symmetries does not imply that phase singularities in the medium do not interact: as mentioned above, the simulations were allowed to converge to a state with stationary phase singularities exhibiting periodic dynamics. In this transient time, spiral waves can have interacted with each other or pushed each other out of the medium. All the above analysis means is that the phase singularities have settled down outside their mutual interaction ranges, probably because for this specific model and parameters B1, the spiral waves are repelling each other. In contrast, the formation of pairs or clusters of spiral waves that settle down in a minimum of an attracting potential could have been detected, provided the force survives the linearization of the original equations.

One classical use case for Lyapunov exponents is bifurcation analysis: Bifurcations are qualitative changes of the attractor structure along a certain path in parameter space. Locally in phase space, certain directions change their stability properties from stable to unstable or vice versa. Therefore, exactly at the bifurcation, there are typically one or more additional neutrally stable directions that can be detected via Lyapunov stability analysis (and yield additional zero Lyapunov exponents). A typical and well-understood bifurcation in excitable media is the transition from a rigidly rotating spiral to a meandering spiral, which has been shown to occur via a Hopf bifurcation [13–15] introducing a second frequency into the system. D. Barkley obtained this important finding by carrying out a linear stability analysis on the fixed point representing a rigidly rotating spiral in the co-rotating frame of constant angular velocity. In this frame of reference, two eigenvalues of the Jacobian (at the fixed point) cross the imaginary axis at the critical bifurcation parameter value.

Figure 3.9 shows the Lyapunov spectrum for a single spiral wave in the Barkley model using parameter set B3. The transition to meandering occurs from right to left for a parameter value of and is identical to the one considered in a publication by Barkley [14]. Towards the bifurcation point (from the right), two Lyapunov exponents approach zero, corresponding to the two-dimensional eigenspace for the fixed-point problem considered by Barkley. At the bifurcation, the rigid rotation of the spiral becomes unstable via the above-mentioned Hopf bifurcation. Since Lyapunov spectra only characterize the attractor the system is currently exploring, the two Lyapunov exponents do not become positive (as the real parts of the eigenvalues for the fixed-point analysis do). Instead, after the bifurcation (towards the left), there are now four zero Lyapunov exponents: two for the periodic or quasi-periodic dynamics on a torus and two—as before—as a result of the Euclidean symmetry (see Sects. 3.1.2 and 3.1.3). The Lyapunov vectors belonging to the two exponents approaching zero towards the critical parameter value correspond to the bifurcating directions in tangent space. Thus, the two vectors shown in Fig. 3.10 span the two-dimensional bifurcating subspace. As the Lyapunov vectors depicted in Fig. 3.4, which represent the Euclidian symmetry, the Lyapunov vectors for the meandering modes are concentrated at the wave front and back. However, while the symmetry modes have an angular periodicity of the wave front modulation of approximately (because of its action as a spatial shift), the meandering modes exhibit a higher spatial frequency. Note that this is similar to the increasing spatial frequency for higher-index modes for plane waves as depicted in Fig. 3.1.

Of course, Lyapunov stability analysis cannot replace linear stability analysis of the corresponding fixed point problem: the most important disadvantage is that the imaginary parts of the eigenvalues cannot be reconstructed from the Lyapunov exponents. Furthermore, as mentioned above, Lyapunov analysis is always restricted to the attractor the system is currently evolving on, while a linear stability analysis can determine the stability properties of arbitrary points in phase space, no matter whether they are stable or not.³ The latter might, however, also be considered an advantage, since the fixed points (or other manifolds) do not have to be determined separately from the stability analysis. As an example, the angular wave speed of the co-rotating frame along with the fields of , and of the Fenton-Karma model are obtained using a Newton-Raphson method in Ref. [16], and the dominant eigenvalues and eigenmodes at the fixed point have to be determined by further sophisticated methods. In contrast, Lyapunov exponents and vectors can be determined by standard procedures without any further algorithmic effort. Note also that the ability to detect and characterize this particular meandering bifurcation via linear stability analysis is based on the geometrically simple trajectory from which it originates, which makes it easily possible to transform a periodic orbit into a fixed point. While this can in the general case always be done by defining appropriate Poincaré sections, switching from the continuous to a time-discrete system might hide important details. In contrast, Lyapunov analysis could also be applied for meandering transitions where the (non-meandering) spiral trajectory is not just a circular motion.

So far, the stability analysis has been applied to homogeneous excitable media, only. However, an important problem in the context of the control and development of undesired activation patterns in excitable media is how heterogeneity influences their stability (see introduction, Sect. 1.​6). One more specific question one might ask is: If a spiral wave is prone to some sort of instability, does heterogeneity stabilize or destabilize it?

The scenario in which this question will be answered uses a specific parameter range in the Barkley model, in which and are fixed, whereas the excitation threshold is varied. Figure 3.11 shows the Lyapunov spectra for a spiral wave in a homogeneous medium with the different values of . For the whole range, the spiral wave rotates steadily, indicated by three Lyapunov exponents which are equal to zero (similar to the spiral in Sect. 3.1.4 before the meandering bifurcation, see Fig. 3.9). The two Lyapunov exponents associated with meandering come very close to zero towards , but instead of bifurcating, they stabilize again for further increasing . Finally, for very high , the medium loses its ability to form spiral waves and so the Lyapunov spectrum becomes all negative, indicating a stable fixed point which corresponds to the stable, homogeneous resting state of the medium.

The heterogeneity that will be considered is a circular region of radius close to the spiral tip, in which the excitation threshold is increased from (which is used in the rest of the medium) to . It is well established that spiral waves are usually drawn to regions of longer spiral wave periods or lower excitability [17–19], although exceptions are known [20]. For the parameters used here, the expected behavior is indeed observed and the spiral wave is attracted by the heterogeneity if initiated close to it. A snapshot of the long time behavior for an exemplary simulation is shown in Fig. 3.12a. From Sect. 3.1.2 it is known that usually, two Lyapunov exponents which are equal to zero indicate the translational symmetry of the activation pattern. Consequently, for a heterogeneity which pins the spiral wave to a specific position as in Fig. 3.12a, these two exponents should become negative, similar to the pinning effect of the numerical grid observed in Fig. 3.5.

However, a parameter scan for and reveals that the behavior of the system is much richer than expected. Figure 3.12c shows the Lyapunov spectrum of the dynamics for a small heterogeneity of in model length units, plotted against the parameter within the heterogeneity. corresponds to a homogeneous medium. As expected, the two Lyapunov exponents and indicating spatial symmetry become negative with increasing excitation threshold . The heterogeneity thus seems to create one specific favorable location for the spiral and the “pinning force” to this location increases with the excitation threshold (indicated by the more and more stable symmetry modes).⁴ One can already see in the same plot that the destabilization of the meandering modes which was visible for homogeneously increased (see Fig. 3.11) is not fully suppressed: and , which were identified as the meandering modes in Sect. 3.1.2, are slightly destabilized, although not enough to lead to a bifurcation and less than in the homogeneous case (see Fig. 3.11). However, if the heterogeneity is larger, as is the case for Fig. 3.12d (), exactly this happens. Not only is the rate with which the symmetry modes become stable much faster, but at the same time, the meandering modes are destabilized much stronger, leading to a bifurcation close to which results in one additional zero Lyapunov exponent above (similar to meandering in a homogeneous medium). Note that, for this to happen, the Lyapunov exponents which correspond to the symmetry and meandering modes have to reverse their order (see value of marked by “” in Fig. 3.12d). This will be discussed in detail at the end of this section. Visual inspection of the spiral tip trajectories confirms that up to , the spiral wave is indeed held in a fixed position by the heterogeneity, as seen in Fig. 3.12a, whereas above the bifurcation value, the tip starts to meander. This effect is seen more easily for an even larger heterogeneity (), with a Lyapunov spectrum as plotted in Fig. 3.12e. The corresponding spiral tip trajectory for the meandering region with two zero Lyapunov exponents is shown in Fig. 3.12b.

For small , the heterogeneity therefore acts as a pinning center, despite the fact that a homogeneous medium of the same increased excitation threshold is close to a meandering instability (cf. Fig. 3.11). For larger , however, this potential instability does not only reappear but is even enhanced, such that an actual bifurcation can take place which is not observed in the homogeneous case. Another interesting phenomenon for is that the symmetry modes are not stabilized for moderate values of , indicating that within the heterogeneity, the spiral can still evolve freely if the heterogeneity is not too strong. Also, note that the Lyapunov spectra for moderate (below the bifurcation value) look very similar to the homogeneous case in Fig. 3.11. The value of at which the spiral wave feels that the modification of is not global and deviates from the homogeneous behavior therefore increases with increasing . This is consistent with the fact that the case of is identical to the homogeneous scenario of Fig. 3.11.

For a stable spiral wave in a homogeneous medium it was shown in Sect. 3.1.2 that the leading Lyapunov spectrum consists of three neutrally stable modes which represent a perturbation along the trajectory of the system and spatial shifts in the two available directions due to Euclidian symmetry. The fourth and fifth Lyapunov exponent were associated with the modes which were found to be responsible for the meandering bifurcation considered in Sect. 3.1.4. However, as seen in this section, the relative order of the symmetry and meandering modes in terms of growth rate is not universal, but can be reversed by interaction of the spiral wave with a heterogeneity. In Fig. 3.12d, it is seen that this reversal can take place when both modes are actually stable, which makes it impossible to differentiate them by observing the dynamics of the spiral wave. By following the change of Lyapunov exponents from the homogeneous case () to the value of indicated by “”, it was implicitly assumed that the trend in growth rate for both modes continues and that indeed the Lyapunov exponents of the meandering modes approach a value of zero (later confirmed when the bifurcation takes place at about ). However, since the ordering of Lyapunov exponents from largest to smallest alone defines their indices, one would not be able to follow this (heuristic) reasoning if only the spectrum at (“”) was known. Fortunately, this is possible making use of one quantity disregarded so far this section: the Lyapunov vectors. Figure 3.13 shows the Lyapunov vectors for the five leading exponents in the system with and . The first one corresponds to the perturbation along the trajectory, whereas the spatial symmetry modes can be found as the fourth and fifth Lyapunov vector. represents a shift of the wave pattern towards the lower right, whereas shifts the spiral wave towards the lower left. In contrast, the second and third Lyapunov vector exhibit a higher-frequency modulation of the wave front characteristic for the meandering modes (see Sect. 3.1.4). In this way, the reversal of stability order necessary for the pinned meandering bifurcation observed in this section can be objectively confirmed.

One of the advantages of Lyapunov stability analysis compared to the linear stability analysis of fixed points is that it can be applied to arbitrary trajectories. Therefore, the heterogeneity in the system does not have to be “simple” for the method to work. Thinking about the complex system of the heart and its inherent heterogeneity, not only abrupt changes of tissue properties in localized regions are conceivable, as considered in Sect. 3.1.5. Rather, real-world (and in particular biological) systems are usually subject to natural variability which leads to a random distribution of parameters. Following this line of thoughts, the excitation threshold of the medium is now modified randomly. To compare the results with those in Sect. 3.1.5, random variations of from the base value of towards higher excitation thresholds are imposed, where a well-defined correlation length serves as a similar measure as the radius of the heterogeneity in Sect. 3.1.5. To construct such a randomly varying , so-called Gaussian random fields are utilized. These are fully defined by their auto-correlation function, which makes it particularly easy to obtain different realizations. If is the intended correlation length, let
be the imposed (spatial) auto-correlation function. Via the Wiener-Khinchin theorem, the Fourier transform of is the power spectrum of the underlying field , so
where denotes the two-dimensional Fourier transform in space. As the Fourier transform of the field is only determined up to its amplitude by this relation, different realizations of the field can be obtained by choosing random phases for each mode (numerically: for each computational point) and transforming the spectrum back to real space:
where is the real part. The standard deviation of the obtained realization across the whole domain is then normalized to one and negative regions are set to zero to give the final modulation :