Terminating transient chaos in spatially extended systems Open Access

Chaotic dynamics occurs in many theoretical, experimental, and real-life systems, and the underlying dynamics has been studied extensively in the past few decades. A particular incentive for this field of research is to control the dynamics in such systems.¹ Various concepts and techniques have been developed for situations when the emerging chaotic dynamics is undesirable. Major breakthroughs were the development of the Ott, Grebogi, and Yorke (OGY) control^2,3 and the Pyragas control,⁴ which both enabled the stabilization of unstable period orbits embedded in the chaotic attractor. Besides stabilizing the dynamics, the question of how to reach a specific neighborhood inside the attractor, called “targeting,” was also investigated in the 1990s.^5,6

However, in most studies dealing with the control of chaotic systems, the underlying chaotic dynamics was persistent, implying the presence of a chaotic attractor. Only some studies consider systems exhibiting transient chaos, where the dynamics is mainly governed by a chaotic saddle. Transient chaos occurs in diverse systems⁷ as, for instance, in ecology,⁸ turbulence,^9–11 coupled semiconductor oscillators,¹² neural networks,^13–15 or NMR-lasers.¹⁶ Chimera states may also often be considered as (long living) chaotic transients,¹⁷ and their collapse has recently been related to epileptic seizures.¹⁸ Also, the chaotic dynamics of spiral/scroll wave dynamics underlying cardiac arrhythmias can be transient.^19,20

Since unstable periodic orbits are also present in the case of transient chaos, the OGY control and Pyragas control methods can be applied too²¹ (although a real chaotic attractor does not exist). Stabilizing the chaotic dynamics in the transient case has been considered, for example, in the case of the tent map,²² the Hénon map,²³ electrical power systems, ecology, or chemical reactions.²⁴ 

In other cases of transient chaos, the chaotic behavior is desired or even essential for the function of a system and the transition toward a stable attractor needs to be prevented. This concept called “chaos maintenance” is relevant, for instance, in the case of (unperturbed) species extinction²⁵ and has been studied in the case of the Lorenz system too.²⁶ 

However, the concepts and methods discussed so far have in common that the state of the system remains inside the chaotic attractor (persistent chaos) or the state space region determined by the chaotic saddle (transient chaos). Let us now assume that these regions of the state space are in general not beneficial and thus it is highly desirable to transfer the system toward another attractor. In the case of persistent chaos, the task of terminating the chaotic dynamics is often discussed and investigated in the context of “multistability”²⁷ and “attractor hopping.” However, this approach is only applicable if the dynamics is governed by an actual attractor, which is not the case with transient chaos. In Sec. I A, we will discuss why this distinction between persistent and transient chaos has a crucial impact on how to terminate the chaotic dynamics of such systems.

In many situations, the chaotic dynamics observed in a system is not desired like the so called “shimmy wheel” in engineering.²⁸ In fact, in certain cases, it can be life-threatening as during cardiac arrhythmias like ventricular fibrillation. As long as the chaotic dynamics is monitored using any (easily accessible) observable, one cannot (in general) distinguish whether the dynamics is persistent or transient. That is why distinguishing persistent chaos from very long transients is technically difficult. However, the state space structure underlying the dynamics has very different properties (chaotic attractor vs chaotic saddle). If the chaotic dynamics is not desired, a perturbation can be applied with the aim of kicking the system out of the chaotic attractor and its basin of attraction. In this context of applying finite perturbations, the concept of basin stability plays a crucial role.^29–33 

Due to the spatial extent of the attractor and its basin [illustrated in Fig. 1(a)], perturbations of a finite size might be insufficient to exit the basin of attraction, no matter in which direction in the state space they are applied [lower yellow arrow in Fig. 1(a)]. Therefore, perturbations of a certain size only can kick the system out of the basin of attraction and let the system evolve to the desired state. In the transient case [Fig. 1(b)], almost all states will converge to the final (desired) state (in the simplified case of only one additional attractor) after a specific transient time (which can, however, be unacceptably long). Due to the chaotic nature of the dynamics, the transient times of states, which are in the vicinity of the reference point in state space, are usually broadly distributed.^34,35 Hence, there are nearby points with larger transient times, but also, most probably neighboring points associated with a significantly lower transient time than the one corresponding to the current state. This is true, in particular, for high-dimensional systems. By perturbing the current state to such a neighboring point in state space [yellow arrow in Fig. 1(b)], the transient time can be reduced significantly and the end of the chaotic dynamics is achieved much faster.

Eventually, the fundamental difference between the termination of chaotic dynamics in the case of persistent chaos (chaotic attractor) and transient chaos (chaotic saddle) is that in the latter case, there is in general no lower boundary for the size of the perturbation, which is necessary to kick the system toward a trajectory quickly converging to the desired state. Hence, although the difference between transient and persistent chaos may not be apparent during the chaotic part of the dynamics, it can play a major role when termination is desired.

In Sec. I A, we explained the general idea of reducing the transient time by the application of a proper perturbation. Finding a proper perturbation in high-dimensional (spatially extended) systems, however, is a costly and elaborate task. Therefore, we use a control scheme originally developed by Cornelius et al.³⁶ for the control of networks. In this study, we apply the control algorithm to two spatially extended systems, which exhibit transient chaos, the Fenton–Karma model³⁷ describing the propagation of electrical excitation waves inside the heart muscle and the Morris–Lecar neuron model³⁸ describing the ion currents of the barnacle giant muscle fiber.

In this section, we will describe the models we consider in this study and explain the control algorithm in a general form afterward.

In numerical simulations of the Fenton–Karma model³⁷ in a two-dimensional rectangular domain (⁠$Lx×Ly=20mm×20mm$⁠), we investigate chaotic spiral waves dynamics [Fig. 1(a)]. The spatial domain was discretized by $Nx×Ny=100×100$ grid points. The dynamics is determined by a system of reaction–diffusion equations (1) and (2), where the dynamics of the secondary variables $h=(v,w)$ and the exact form of the currents $Ifi$⁠, $Iso$ and $Isi$ are shown in the supplementary material,

The values of all parameters including the diffusion constant $D$ and the membrane capacitance $Cm$ [in Eq. (1)] are also given in the supplementary material. An exemplary snapshot of the spiral wave dynamics is shown in Fig. 2(a).

For the second model, we consider a one-dimensional ring of $N=100$ Morris–Lecar neurons. The dynamics of the membrane potential of a neuron is determined by Eq. (3), where $D$ is the diffusion constant, $Cm$ the membrane capacitance, and $IL$⁠, $ICa$⁠, and $IK$ denote transmembrane currents, respectively. The behavior of the number of open potassium channels $ni$ is given by Eq. (4). The form of the transmembrane currents, $nss(Vi)$⁠, $τn(Vi)$ and simulation parameters are given in the supplementary material,

The algorithm used for the control of the models discussed above is based on an approach, which was developed by Cornelius et al.³⁶ for the control of networks. The ultimate task of the algorithm is the following: assuming that the system is in the current state $xstart(t0)$ at $t0$ , we look for a perturbation $ptotal$ to the state, such that the perturbed state $xstart$ + $ptotal$ evolves within a fixed period $Tevo$ to (or close to) the desired state $xdesired$⁠,

We parameterize $ptotal$ such that for both models, it acts only on the first variable of the system [$Vm$ (Fenton–Karma) and $Vi$ (Morris–Lecar), respectively]. $ptotal$ consists of a combination of individual spatially localized perturbations $pi$⁠. In the case of the Morris–Lecar network, such a single perturbation could be applied to each neuron (⁠$i∈{1,…,100}$⁠), in the case of the Fenton–Karma model (simulation domain of $Nx×Ny=100×100$ grid points) single perturbations had a size of $2×2$ grid points (resulting in $i∈{1,…,2500}$ in order to cover the whole simulation domain). The perturbation amplitude of a single perturbation is set to $ΔVm=0.5a.u.$ (Fenton–Karma) and $ΔVi=0.1mV$ (Morris–Lecar), respectively. Despite this parameterization, the number of possible perturbations $ptotal$ is due to the high (spatial) dimensionality of the systems still considerably large. Therefore, the main purpose of the algorithm is to suggest which of the single (local) perturbations $pi$ should be included into $ptotal$⁠.

At the beginning, a value for the period $Tevo$ has to be chosen, which specifies a time period until when the dynamics should reach the desired state. The initial state is equal to the current state of the system: $x0=xstart(t0)$⁠.

1. Apply each single perturbation to the initial state individually: $xperti=x0+pi$ and evolve each perturbed state in time for $Tevo$ [see Fig. 3(a)].

2. Out of all possible single perturbations $pi$⁠, choose the perturbation whose trajectory comes closest to the desired state (based on a metric $d(x,y)$⁠):

   $pclosest=minpi,t∈[0,Tevo]d(xperti(t),xdesired).$

   (6)

3. Add this chosen single perturbation renamed as $pcj$ (where $j$ is the number of iterations of the algorithm) to $ptotal$⁠.

4. Update the initial state $x0→x0+pcj$⁠.

5. Repeat steps 1–4, until in step 2, the metric indicates that by adding the single perturbation $pclosest$ the system reaches the desired state [see Fig. 3(b)]. By finally adding the latest single perturbation $pcj$ (or the current $pclosest$⁠, respectively) to $ptotal$⁠, the final form of $ptotal$ is

   $ptotal=pc1+⋯+pck,$

   (7)

   where k is the number of iterations and thus the total number of single perturbations. $ptotal$ fulfills then the desired property equation (5). Note that in step 1, only perturbations that are not included in $ptotal$ so far are considered (otherwise perturbations would be applied twice at the same spot).

The choice of the metric $d(x,y)$ is crucial for the algorithm, since it is mainly responsible for the choice of single perturbations and thus how to most efficiently navigate through the high-dimensional state space. Hence, we believe that for a specific system, the metric should be adapted regarding the actual dynamics of the system. Using the Euclidean distance as a metric, for example, could in spatially extended systems not necessarily be the most effective choice (measured in terms of an average low number of single perturbations). For instance, the main governing objects of the dynamics of Fenton–Karma simulations are a relatively low (compared to the actual dimension of the system $100∗100∗3$⁠) number of spiral waves. The ultimate aim of reducing the number of these spiral waves (computed by the number of phase singularities) is somewhat compatible with bringing each cell toward the resting state. Both strategies would have the same (or nearby) desired point(s) in state space, still the individual trajectories resulting from differently chosen single perturbations would probably significantly differ. In the case of Fenton–Karma simulations, we first chose in step 2 of the algorithm those single perturbations which lead to the lowest number of phase singularities (and thus spiral waves). If this choice was not unique, we decided to keep the single perturbation $qi$⁠, which additionally had the lowest Euclidean distance to the resting state. In the case of the Morris–Lecar network, however, due to the lack of a better alternative, the Euclidean distance was used for $d(x,y)$⁠.

We applied the control algorithm with varying $Tevo$ for 100 initial conditions for each model. In order to exclude that due to the transient nature the unperturbed trajectories would self-terminate anyway, we ensured that the self-termination does not occur for at least 20 s (Fenton–Karma) and 50 s (Morris–Lecar), respectively. After demonstrating a working example, we investigate the dependence on the parameter $Tevo$⁠, compare the control scheme in terms of efficiency to random perturbations and discuss possible underlying mechanisms of the control scheme. The results shown here contain simulations of the Fenton–Karma model only. The associated analysis regarding the Morris–Lecar network provided very similar results compared to the Fenton–Karma model, which is given in the supplementary material.

An exemplary course of the algorithm for an episode of spiral wave chaos is depicted in Fig. 4 (Multimedia view).

We followed the control scheme with a chosen evolution period of $Tevo=500ms$⁠. The scheme finished after three iterations and provides the perturbation $ptotal=pc1+pc2+pc3$⁠. The perturbed trajectory [starting with $xstart+ptotal$ depicted in (a) at $t1=0ms$] terminates within $Tevo$⁠, whereas the spiral wave dynamics of the unperturbed trajectory [subplot (b)] continues. The difference between the perturbed and unperturbed trajectory $δVm$ [shown in subplots (c)] demonstrates how the initially spatially localized single perturbations $pci$ affect a growing part of the whole simulation domain as time evolves. Note that in the first plot of Fig. 4(c) (Multimedia view), the pure perturbation $ptotal$ is depicted, consisting in this example of three spatially localized perturbations. For the sake of clarity, we remind the reader that different choices of $ptotal$ would result in different combinations of these spatially localized perturbations.

The initial choice of the temporal evolution period $Tevo$ is essential for the control scheme. It turns out that the larger we choose the length of $Tevo$⁠, the less perturbations are needed in order to control the dynamics. We investigated 100 initial conditions and applied the control scheme with varying $Tevo$⁠. The resulting statistics of how many of the initial conditions could be controlled using $Npert$ single perturbations is shown in Fig. 5(a). Decreasing the length of $Tevo$ shifts the distribution of successfully controlled cases toward an increased $Npert$ and increases also the number of initial conditions, which require more than 10 single perturbations. The accumulated amount of initial conditions, which could be controlled with less or equal than 10 single perturbations, is shown in Fig. 5(b). Note that with $Tevo≥800ms$⁠, all 100 initial conditions could be controlled with $Npert≤10$ single perturbations.

The unperturbed dynamics is transient itself, meaning each trajectory will reach the resting state at some point without any interaction. So, just by randomly perturbing the current trajectory, we also change the transient time. That is why we demonstrate that the control scheme is nonetheless highly efficient, in terms of finding “shortcuts” in the neighborhood of the current state that lead as quick as possible to the desired state. For this purpose, a fixed number of single perturbations (shape and amplitude are the same as before) were distributed randomly (10 different distributions) and were then applied to the 100 initial conditions that were used before. Thus, per fixed number of perturbations, $10∗100=1000$ simulations were performed and the success rate was computed (in terms of whether the system reached the resting state within 500 ms). In this way, a termination rate could be computed, dependent on the number of randomly distributed single perturbations. The resulting curve is shaped sigmoid-like in Fig. 5(c). It is noteworthy that at a number of $Npert=500$ single perturbations, the termination rate is below 7%, whereas an exemplary pattern corresponding to $Npert=500$ single perturbations is depicted in Fig. 5(d). Comparing the differences in magnitude of order of the x-axes of Figs. 5(a) and 5(c), respectively, and the significant difference between the perturbation patterns in Figs. 5(d) and 4(a) (Multimedia view) demonstrates the efficiency of the control scheme used above. In the case of randomly distributed perturbations, a success rate of 50% is achieved only at approximately $Npert=842$⁠.

In Secs. III A, III B, and III C, we demonstrated the operability and performance of the control scheme. However, although we showed that the single perturbations suggested by the algorithm are highly efficient in finding “shortcuts” toward the desired state, the specific distribution of the perturbations was determined with a brute force-like method. In terms of general control theory, it is of great interest to understand these distributions. As a first approach, we studied the impact of single perturbations on the system one by one. For an exemplary trajectory of the Fenton–Karma model, the unperturbed initial state $xstart$ is shown in Fig. 6(a), together with consecutively added single perturbations, which were iteratively determined by the control scheme.

In subplot Fig. 6(b), each state is evolved in time for $Tevo=500ms$⁠. By adding single perturbations to the initial state, the final state (lower row) changes in parts of the domain. However, single perturbations seem to act on specific regions of the domain, only. For example, by adding the second perturbation $pc2$⁠, the spiral wave in the lower part (marked by the red circle) has vanished. The rest of the spatial domain, however, seems to be unaffected by the second perturbation. Only adding the third perturbation $pc3$ allows for the entire system to reach the desired resting state. The observation that single perturbations mainly affect spatial domains, which are confined to the neighborhood of the perturbation’s origin, is closely related to the diffusive nature of the dynamics. From this point of view, a possible approach to understand the efficient control of these systems could imply that the spatial domain needs to be divided into meaningful segments, and each segment needs to be controlled individually. In fact, in the case of spiral wave dynamics in excitable media, Byrne et al. suggested to divide the spatial domain into “tiles” based on the governing dynamics of spiral waves.³⁹ 

We demonstrated and investigated in this study the fundamental difference between the termination of persistent and transient chaos, respectively. In the first case, a perturbation to the current state of the system is applied with the aim to escape from the attractor’s basin. In the latter, the aim of the perturbation is to significantly reduce the transient time by finding a nearby state that converges to the desired state much faster. This difference leads to the conclusion that in the presence of a chaotic attractor, a perturbation of a minimum size is required, in order to kick the system out of the attractor and its basin. For transient chaos, there is in principle no such lower boundary. Thus, although the difference between transient and persistent chaos cannot be perceived by observing the chaotic dynamics (unless self-termination has occurred), it can play a major role when the termination of the chaotic dynamics is required.

We investigated two spatially extended systems that exhibit transient chaos. The Fenton–Karma model describes action potential propagation in cardiac tissue, and the Morris–Lecar network is a model for neuronal activity. Using an algorithm originally developed by Cornelius et al.,³⁶ we were able to terminate the dynamics using a small number of spatially localized perturbations.

Terminating chaotic dynamics with a small perturbation only is relevant also in real-life systems like the defibrillation of life-threatening cardiac arrhythmias, whose dynamics is governed by chaotic spiral/scroll wave dynamics.^40–43 Since side effects of conventional defibrillation increase with the energy of the electrical shock (corresponding to the size of the perturbation), the aim is to significantly reduce the amount of energy.^44–46 Although great progress has been made in the past by reducing the energy and maintaining a high success rate,^44,47 the question whether there is a natural lower energy boundary for successful defibrillation is not answered yet. The existence of such a minimal energy that is inevitably necessary to successfully terminate the chaotic spiral/scroll wave dynamics would indicate the presence of a chaotic attractor with a basin of a certain size [Fig. 1(a)]: a minimum defibrillation energy (perturbation size) is required in order to escape from the attractor and its basin. However, since there are indications that the underlying spiral/scroll wave dynamics can be transient^19,20,48 and cases of self-terminating arrhythmias have been reported,^49–51 the above-mentioned considerations suggest that such a lower boundary does not exist (at least for a subset of arrhythmias). However, in order to more efficiently terminate such a dynamics, more sophisticated and complex perturbation patterns/techniques are required. Furthermore, the efficient termination of transient chaos can also be desired in other (real-life) spatially extended systems, such as turbulent flow,^10,11 or in chemical reactions.⁵² 

The control algorithm applied can so far not be easily transferred to real-life (experimental) systems. Nevertheless, this study demonstrates as a proof of principle that control and targeting of high-dimensional systems exhibiting transient chaos can be achieved with exceptionally small interactions with the system and will hopefully stimulate further research on coping with transient chaos in real world applications.

See the supplementary material for additional technical details and the results regarding the Morris–Lecar network.

We acknowledge the International Max Planck Research School (IMPRS) for Physics of Biological and Complex Systems (PBCS), the Federal Ministry of Education and Research (BMBF, No. FKZ031A147, GO-Bio), and the German Centre for Cardiovascular Research (DZHK) e.V. for financial support. Furthermore, we thank Stefan Luther for fruitful discussions and continuous support.

The data that support the findings of this study are available from the corresponding author upon reasonable request.