Spontaneous termination of chaotic spiral wave dynamics in human cardiac ion channel models

Chaotic spiral or scroll wave dynamics can be found in diverse systems. In cardiac dynamics, spiral or scroll waves of electrical excitation determine the dynamics during life-threatening arrhythmias like ventricular fibrillation. In numerical studies it was found that chaotic episodes of spiral and scroll waves can be transient, thus they terminate spontaneously. We show in this study that this behavior can also be observed using models which describe the ion channel dynamics of human cardiomyocytes (Bueno-Orovio-Cherry-Fenton model and the Ten Tusscher-Noble-Noble-Panfilov model). For both models we find that the average lifetime of the chaotic transients grows exponentially with the system size. With this behavior, we classify the systems into the group of type-II supertransients. We observe a significant difference of the breakup behavior between the models, which results in a distinct dynamics during the final phase just before the termination. The observation of a (temporally) stable single-spiral state affects the prevailing description of the dynamics of type-II supertransients as being “quasi-stationary” and also the feasibility of predicting the spontaneous termination of the spiral wave dynamics. In the long term, the relation between the breakup behavior of spiral waves and properties of chaotic transients like predictability or average transient lifetime may contribute to an improved understanding and classification of cardiac arrhythmias.


Introduction
Transient chaos is a widespread phenomenon, where the chaotic dynamics of a system is not persistent but decays after some time. Chaotic transients appear in diverse systems [1] as, for instance, in ecology [2], turbulence [3], coupled semiconductor oscillators [4], neural networks  [5][6][7] or NMR-lasers [8]. Regarding the control of such chaotic transients it depends on the specific system or objective whether the chaotic phase should be extended or a quick self-termination is desired. An example for the latter case can be found in the field of cardiac dynamics, where chaotic spiral or scroll waves determine the electrical excitation dynamics during life-threatening cardiac arrhythmias like ventricular fibrillation [9][10][11][12]. Spontaneous termination of such ventricular fibrillation in human patients has been reported [13] and it has been discussed how the underlying mechanisms of self-terminating arrhythmias could improve the understanding and efficiency of antiarrhtyhmic drugs [14,15]. Also, in numerical simulations it was found that chaotic spiral wave dynamics can in fact be transient [16][17][18]. However, it was not yet clear whether the results hold also for human cardiac ion channel models. Furthermore, former studies showed that chaotic transients of spiral and scroll wave dynamics observed in excitable media can be assigned to the group of type-II supertransients, which are characterized by a spontaneous self-termination, which is, so far, not predictable a significant amount of time ahead. Although it has been shown that the state space structure itself changes nearby the "exits" of the chaotic regime [19,20], precursors based on conventional variables have not been found yet.
In this study we show that the transient nature of chaotic spiral wave dynamics is a robust phenomenon which holds also for human cardiac ion channel models, explicitly the Bueno-Orovio-Cherry-Fenton model [21] and the Ten Tusscher-Noble-Noble-Panfilov model [22,23]. Furthermore, we reveal that the wave breakup behavior of the dynamics can have a large impact on the mechanism that leads to self-termination, and in this way may have implications for a possible prediction of self-termination of spiral wave dynamics, which we observe during ventricular fibrillation.
This study is organized in the following way: In the next section we present the investigated models and describe how we determined average transient lifetimes hTi in those systems. In the next section we demonstrate that the average lifetimes increase exponentially with the system size in both models, whereas the number of spiral waves increases linearly. Then, we investigate how the differences of spiral wave breakup behavior between both models can affect the final episode of the transients. In the last section, we discuss how these findings can contribute towards an improved understanding of the underlying mechanisms of spiral wave chaos, which occurs, for instance, during ventricular fibrillation and how future simulations of cardiac arrhythmias could benefit from these results.

Materials and methods
We investigate two models of cardiac ion channel dynamics in a two-dimensional squared shaped domain with the domain size A = L x × L y , where L i is the length of the domain in direction "i". The underlying reaction-diffusion equations of both models will be explained in the next section. The creation of initial conditions and the detection of phase singularities is discussed in the following section. Further details about both methods and a detailed description of the algorithms are given in S1 Appendix.

The system of reaction-diffusion equations
The propagation of electrical excitation waves through cardiac tissue was modeled by a system of reaction-diffusion equations describing the evolution of the transmembrane potential V m : where C m is the membrane capacitance, D the diffusion tensor (in our case isotropic), I ion the net transmembrane ion current, and I stim an external stimulation current, which is used for the creation of initial conditions only (details about the creation of initial conditions are given later). The individual transmembrane currents within I ion as well as additional ion channel dynamics constitute the remaining body of differential equations. In this study we investigated the Bueno-Orovio-Cherry-Fenton model (from now on denoted as the BOCF model) and the Ten Tusscher-Noble-Noble-Panfilov model [22] with the update of the original model regarding the description of the calcium dynamics, as discussed in [23] (from now on denoted as the TNNP model).
The Laplacian (first term on the right-hand side of Eq (1)) was discretized using a ninepoint stencil. The temporal integration was performed through a first-order explicit Euler step for all equations except those pertaining to the dynamics of gating variables, which were evolved through a first-order Rush-Larsen scheme instead [24]. Most computations were carried out on a cluster of 32 Nvidia Tesla P100 GPUs with 16 GB VRAM and a 32-bit floatingpoint unit performance of 10 TFLOPS each.
The Bueno-Orovio-Cherry-Fenton model. Instead of modeling the large variety of distinct ion channel currents present in a cardiac cardiomyocyte, the Bueno-Orovio-Cherry-Fenton model aims at describing the net transmembrane current by three contributions: a fast inward current I fi , slow inward current I si , and slow outward current I so . Due to this "effective" description of the dynamics, it is convenient to introduce the dimensionless membrane potential u = (V m − V rest )/(V fi − V rest ) which is normalized by the resting membrane potential V rest and some upper limit, e.g. the Nernst potential of the fast inward current V fi . Hence similarly rescaled currents can be introduced by, for instance, J fi = I fi /(C m (V fi − V rest )). Using this rescaling Eq (1) becomes with the three "effective" currents J fi , J so , and J si : where the rescaled currents have the unit of inverse time. J stim indicates the (rescaled) external stimulation current, which is used for the creation of initial conditions. In the Bueno-Orovio-Cherry-Fenton model (BOCF model) the ionic currents are described by where v, w, and s are gating variables, θ v , u u , u o , and θ w model parameters, and τ fi , τ o , τ so , and τ si time constants. A value of D = 2 × 10 2 mm 2 /s was selected for the isotropic diffusion tensor. The model parameters chosen correspond to the "TNNP" parameter set in [21], thus the parameter set aims at reproducing the dynamics of the original TNNP model [22]. For the numerical integration scheme, a spatial resolution of h = 1 mm along with a time step of length Δt = 1 × 10 −4 s was used. A snapshot of the normalized membrane potential u during a typical episode of spiral wave chaos is shown in Fig 1(a).
The ten Tusscher-Noble-Noble-Panfilov model. The revision of the Tusscher-Noble-Noble-Panfilov model (TNNP model) was used for all computations [23]. It features twelve transmembrane currents and six state variables governing the temporal evolution of ion concentrations. The model was implemented based on the description by the CellML project [25] with parameters mostly taken from the epicardial parameter set from [23] with the exception of C m = 0.2 × 10 −2 μF/mm 2 and D = 1.54 × 10 2 mm 2 /s (isotropic). The numerical computations were carried out with a spatial resolution of h = 0.2 mm and a time step of length Δt = 2 × 10 −5 s. Additionally, the "slope 1.8" modifications in [23] were applied to enable spiral wave breakup. A snapshot of the membrane potential V m during a typical episode of spiral wave chaos is shown in Fig 1(b).

Creation of initial conditions and detection of phase singularities
Since relevant quantities discussed in this study vary during a single simulation (e.g. the number of spiral waves) or between different initial conditions (e.g. the transient lifetime), we determined mainly averaged quantities. In order to provide statistically robust statements, it is necessary to create many different initial conditions which sufficiently sample the relevant region of the state space. To create such a set of independent initial conditions, we used a sequence of spatially distributed local stimuli. Details about this process can be found in the Appendix. An exemplary video of the creation of an initial condition of the BOCF model is given as supporting information (S1 Video).
The number of spiral waves was determined via detecting the phase singularities which are present in the system. This procedure was done by a two-step protocol: In a first step, the phase θ was determined by a two-dimensional phase-space projection. Phase singularities (which can be associated with the tips of the spiral waves) could then be identified at a specific position by evaluating closed line integrals of the gradient of θ encompassing the spot in space. Further details about the procedure can be found in the Appendix.

Determination of the average transient lifetime
We quantify chaotic transients by determining their average lifetimes which correspond to a specific domain size A. In a system exhibiting transient chaos, the number of initial conditions with a transient lifetime of at least t from a sample of size N, N Ch (t), is expected to decay exponentially [1]: where κ denotes the escape rate describing the "escape" from the chaotic regime of the state space which is usually determined by a chaotic saddle, towards the non-chaotic region. In our case, the latter can be identified as the global resting state of the system without excitation

Results
The presented results are structured in the following way: In the first section, we discuss how the average lifetime of chaotic spiral wave episodes depends on the system size. In the following part, we show that the average number of spiral waves grows linearly with the system size.
We then demonstrate how the number of spiral waves is distributed differently in both investigated models, respectively, and how this has an impact on the correlation between the average number of spiral waves and the end of the chaotic episode (last section).

Transient lifetime and average spiral wave number depend on the system size
The spiral wave dynamics we observe in simulations of both models is transient, meaning that after a certain amount of time no more spiral waves are present and the system returns globally to the resting state. During the chaotic episode, a fluctuating number of spiral waves can be  (6), where data points below the fit threshold were discarded. We characterize the transient chaotic spiral wave episodes in both models by determining the corresponding average lifetime. For this purpose, 200 initial conditions were created and their average lifetime was determined (details can be found in the methods section). This procedure was realized for varying system size domains . Fig 3(a) and 3(b) show an exponential increase of the average lifetime hTi with the system size L for the BOCF model and the TNNP model, respectively.
Spatially extended systems with transient dynamics whose lifetime grows rapidly with the system size are called "supertransients". More specifically, the group of type-II supertransients are characterized by an exponential scaling behavior of the lifetime of the transients with the system size [26]: In  Table 1.
The spiral waves are the organizing objects which determine the dynamics during chaotic episodes. In particular, the creation and annihilation of spiral waves during the episode is the underlying mechanism for the self-termination in the end. In order to quantify how an increase of the domain size changes the spiral wave dynamics, we furthermore measured the average number of spiral waves for each system size A by determining the organizing centers of the spiral waves, so called phase singularities. Details about the process of determination can be found in the supporting information section.
In Fig 4(a) and 4(b) the average number of spiral waves (averaged over time) hN spiral i t is shown for various domain sizes for the BOCF model and the TNNP model, respectively. It is noteworthy that the subscript t in hN spiral i t indicates that the averaging (denoted by the  (7), which indicates the exponential scaling of the average transient lifetime with the system size. The corresponding fit parameters are given in Table 1. The numerical data underlying this graph and the following figures are also provided as Supporting information S1 Data. Since the number of spiral waves present in the system fluctuates during a typical episode of chaotic spiral wave dynamics, the standard deviation of hN spiral i t indicated by the error bars is rather large. However, a linear increase of the average spiral wave number could be verified in both models (subplots (a) and (b), respectively). This observation is in agreement with similar studies of other ion channel models of transient spiral wave dynamics [16,17].

Wave breakup behavior affects final phase of transients
The underlying mechanism which determines the transient spiral wave dynamics is the continuous creation and annihilation of spiral waves, causing a fluctuating number of spiral waves present in the system, which reaches zero at the point in time of self-termination. Whereas spiral waves can annihilate with each other (pair annihilation) or with the boundary (in the case of no-flux boundary conditions), there are also different processes leading to the creation of spiral waves. C. Marcotte and R. Grigoriev [27] distinguish between "wave breakup" and so called "wave coalescence".
In the first case, the front of an excitation wave interacts with a waveback of the same wave (e.g. induced by a conduction block) and increases the number of excitation waves in this way. "Wave coalescence", however, is denoted as the interaction of an excitation front with a refractory region associated with another wave and the following creation of two separated wavesfronts. Whereas the first mechanism can considered responsible for the transition from a Table 1. Fit parametersã and γ according to Eq (7) linear least-squared fit of the natural logarithms of the average lifetimes hTi depicted in Fig 3(a) and 3(b), for the BOCF model and the TNNP model, respectively). The standard deviation was determined based on the variance of each parameter according to the least-squared method.  single-rotor state (which can be associated with tachycardia) to a more chaotic multi-spiral state (which can be associated with fibrillation), the process of "wave coalescence" is rather linked to the maintenance of an already complex state which exhibits many spirals. In fact, numerical studies have shown that when the system is already in a multi-spiral state, "wave coalescence" is the dominant mechanism for the creation of spiral waves [27] (C. Marcotte and R. Grigoriev did not find a single instance of wave breakup in simulations during multispiral chaos).

Modelã
As also stated in [21] we observe in our simulations that a "wave breakup" occurs only very rarely in the BOCF model, whereas the TNNP model promotes "wave breakup". However, since we initialize both systems in such a way that several spiral waves are present, the number of spiral waves fluctuates in both models. A qualitative difference can only be observed for states of a low number of spiral waves (in particular N spiral = 1). Since wave breakup is not promoted in the BOCF model, the single-spiral state is quite robust and persists for a significant amount of time. This behavior is explicitly relevant during the final phase of the investigated transient episodes. Fig 5(a) and 5(b) show the number of spiral waves N spiral for exemplary episodes of chaotic spiral wave dynamics for the BOCF model and the TNNP model, respectively. Time is normalized such that the respective self-termination occurs at t = 0. In both cases N spiral fluctuates at the beginning of the time series. In Fig 5(a) the system deteriorates into the single-spiral state at around t � −25 s and returns to the resting state much later. However, the transition to the single-spiral state does not necessarily imply a subsequent self-termination. The evolution of N spiral for another episode of the BOCF model is given in S1 Fig, in which the system shows a (temporally) intermediate single-spiral state, but increases its number of spiral waves again afterwards. Exemplary videos of the initial condition discussed in S1 Fig which show the transition from a multi-spiral state into the single-spiral state and the transition back into a multi-spiral state are given as supporting information (S3 and S4 Videos, respectively). In the case of the TNNP model (Fig 5(b)) the general behavior of a pronounced single-spiral state cannot be observed but N spiral fluctuates continuously and reaches zero rather "spontaneously" (see Fig 5(b)).

Quantifying differences of spiral wave dynamics
We want to discuss and interpret the quantitative difference of the spiral wave dynamics between the two models shown in the previous section in a more abstract/mathematical description: As stated in [17], the transient dynamics of spiral wave chaos can be described as a Markov chain, where each state is determined by the number of the spiral waves of the system. Transition probabilities determine how likely a state will increase or decrease the number of spiral waves. For both models, we statistically investigated the probability of finding the system in a state with a specific number of spiral waves, sketched in Fig 6(a) and 6(b), for the BOCF model, and the TNNP model, respectively. For both models we measured N spiral , for three distinct domain sizes for each model, respectively. Except for low numbers of N spiral , a smooth distribution can be observed, with a local maximum which shifts towards larger N spiral for For three different domain sizes each (denoted by different colors) the distribution of the relative appearance of N spiral is shown for the BOCF model (a) and the TNNP model (b), respectively. In both cases, the distribution shifts towards larger N spiral with increasing domain sizes. The shape is, however, biased for low number of spiral waves. In particular in the case of the BOCF model, a distinct peak for the single-spiral state for all three domain sizes can be observed.
https://doi.org/10.1371/journal.pone.0221401.g006 larger domains. However, in the case of the BOCF model the pronounced characteristic of the single-spiral state is clearly visible which is mostly insignificant in the case of the TNNP model (the somewhat non smooth behavior for a domain size of A = 0.66 × 10 5 mm 2 is caused by a relatively small domain size, where the maximum value of the distribution is close to N spiral = 0).
The former investigations show that the wave breakup behavior of the dynamics can have a significant impact on the terminal phase of self-terminating episodes. We verified this observation by measuring N spiral before self-termination for 200 initial conditions for each model, respectively . Fig 7(a) and 7(b) show the number of spiral waves for the BOCF model (domain size of A = 4.10 × 10 5 mm 2 ) and the TNNP model (domain size of A = 1.28 × 10 5 mm 2 ), respectively, averaged over 200 initial conditions hN spiral i IC , where time was normalized such that all initial conditions terminate at t = 0 s. Note that the subscript IC in hN spiral i IC indicates that the average denoted by the brackets is performed over initial conditions, and not over time (in contrast to hN spiral i t shown in Fig 4, for example).
Far from the point of self-termination, hN spiral i IC (thick red line) saturates in both cases at values (hN spiral i IC � 20, and hN spiral i IC � 18, respectively) coinciding with the average number of spiral waves in a system (compare Fig 4(a) and 4(b), respectively). While in the case of the TNNP model (subplot Fig 7(b)) hN spiral i IC starts to deviate significantly from this value around t � À 10s, the value decreases much earlier (around t � À 40s) in the case of the BOCF model (subplot Fig 7(a)). This significant difference between the models is also a consequence of the pronounced characteristic of the single-spiral state before self-termination present in the BOCF model.

Wave breakup behavior of transient spiral wave chaos in the context of type-II supertransients
With the exponential relation between the system size A and the transient lifetime hTi we found the typical scaling behavior of type-II supertransients. In [26] the dynamics of such systems is described as "quasi-stationary" with a spontaneous collapse of the dynamics towards some (low-dimensional) attractor. In the case of chaotic spiral wave chaos, we can typically observe a fluctuating number of spiral waves, and the temporal point of self-termination can not easily be predicted by using "conventional" variables (e.g. N spiral in Fig 5(b)). However, the mechanism of self-termination can in the case of spiral wave chaos be described in an ordered way by using the above-mentioned concept of Markov chains [17]: All trajectories have to pass the state of N spiral = 1 (N spiral = 2) before terminating via collision with the boundary (or via a pair annihilation in the case of N spiral = 2 spirals), described by the Markov chain by a transition towards the state with N spiral = 0. The transitions from one state to another are characterized here by transition probabilities. The governing processes of creation and annihilation of spiral waves allow only a change of the number of spiral waves by ±1 or ±2 (if considering arbitrarily small time intervals). In general, the temporal evolution of a state with N spiral waves from time t to t + Δt can be described by five probabilities: The system can remain in the same state, meaning the number of spiral waves does not change (P ΔN=0 ), or the number of spiral waves can increase, or decrease by different mechanisms (P ΔN=+1 , P ΔN=+2 , or P ΔN=−1 , P ΔN=−2 , respectively). Regarding the Markov chain description, changes of the wave breakup behavior of spiral waves can significantly alter these transition probabilities. In particular for transitions from the N spiral = 1 state, P ΔN=0 is much larger compared to transition probabilities with ΔN 6 ¼ 0 in the BOCF model, as it is in the case of the TNNP model: Since "wave breakup" of a single spiral wave does occur only rarely in the case of the BOCF model, the system remains in the same Spontaneous termination of chaotic spiral wave dynamics in human cardiac ion channel models state for a significant time. For states with N spiral > 1, the number of spiral waves changes much faster (Fig 5(a)) and the transition probabilities with ΔN 6 ¼ 0 are increased (due to wave coalescence).
In the context of type-II supertransients, the single-spiral state before self-termination is a distinct state compared to states with N spiral > 1 (distinguished by significant variations of transition probabilities and occupation times). This does not imply the existence of a precursor of the upcoming collapse in the classical sense, but the temporally structured course of the selftermination (via the "intermediate" N spiral = 1 state) is an observation that extends the "quasistationary" description of the dynamics.

Discussion
We showed in this study that the transient nature of chaotic spiral wave dynamics is not a feature of generic models but is rather a robust phenomenon which can also be found in simulations using human cardiac ion channel models (Bueno-Orovio-Cherry-Fenton model and the Ten Tusscher-Noble-Noble-Panfilov model). The general understanding of cardiac arrhythmias like ventricular fibrillation, whose electrical wave dynamics is determined by spiral or scroll waves, could benefit from the perspective of chaotic transients and its properties: Whereas the number of spiral waves grows linearly with the spatial domain size in both models, the average transient lifetime of the chaotic episodes (which would have a significant relevance when interpreted in terms of the length of cardiac arrhythmias) increases exponentially with the domain size (which might be interpreted as the heart muscle volume). In fact, the role of a critical heart muscle volume which is necessary for sustained arrhythmias has been discussed and also investigated in experimental studies [28,29]. Also, clinical studies show, that the risk for arrhythmias and related mortality and morbidity is increased with a larger heart muscle volume [30][31][32].
Although the scroll wave dynamics during fibrillation of the ventricle is three-dimensional, former studies have shown that the transient features of the spiral wave chaos are also present in the three-dimensional case [16,33]. However, in order to estimate and understand the role of transient spiral/scroll wave dynamics in the context of cardiac arrhythmias, the heterogeneous structure of the cardiac tissue and more realistic heart geometries should be taken into account in future numerical studies.
Novel control methods [34] which aim at reducing the severe side effects of conventional defibrillation methods, could benefit from an enhanced understanding of the governing mechanisms which underlie cardiac arrhythmias and self-termination of its chaotic dynamics.
Furthermore, although chaotic transients of systems characterized by an exponential increase of the average lifetime with the spatial domain size (denoted by type-II supertransients) generally terminate spontaneously and without any obvious indications a significant amount of time before the collapse, we showed that in the specific case of chaotic spiral wave dynamics the wave breakup behavior can have a significant impact on the final phase of an episode. In simulations of the BOCF model we observed that before self-termination the system remains in a relatively stable single-spiral state. This behavior coincides with the clinical observation that recorded self-termination of ventricular fibrillation in human patients showed an intermediate state of higher organization in the electrocardiogram [35,36], e.g. ventricular tachycardia which may be associated with a single rotor. Observing this effect in both, numerical simulations and real hearts, can also be relevant regarding the question what kind of wave breakup behavior predominates during ventricular fibrillation in real hearts (specifically concerning the bidirectional transition between ventricular fibrillation and ventricular tachycardia). This knowledge can be essential when numerical simulations of cardiac tissue should reproduce the electrical wave propagation in a real heart in the most accurate way, and an appropriate models and parameters are required.
Also, in terms of a possible prediction of self-termination, statistically the absolute number of spiral waves present in the system equal to one can in this model serve as an observable indicating a higher chance for an upcoming self-termination.
From the point of view of nonlinear dynamics, the Markov chain description of the dynamics together with the significantly pronounced single-spiral state in the case of the BOCF model extends the prevailing description of type-II supertransient dynamics as "quasi-stationary". This observation might contribute to the challenge of predicting an upcoming self-termination in these systems and the development of related applications in the field of cardiac arrhythmias.
Supporting information S1 Fig. The temporal evolution of the number of spiral waves N spiral during an episode of transient chaos of the BOCF model. During this episode, the system enters the "single-spiral" state (around t � −60 s) but through spiral wave breakup N spiral increases again (around t � −50 s) before it finally self-terminates (at t = 0 s). (EPS) S1 Video. Exemplary episode of the creation of initial conditions. The temporal evolution of the rescaled membrane potential u of a simulation of the BOCF model is shown, during the application of multiple pacing pulses, which were used to create initial conditions.