Suppression of fibrillatory dynamics consisting of stable rotors by periodic pacing

Recent experimental studies have shown that a sequence of low-energy electrical far-field pulses is able to terminate fibrillation with substantially lower per-pulse energy than a single high-energy electric shock (see S Luther et al Nature 475 235–39). During this low-energy antifibrillation pacing (LEAP) procedure only tissue near sufficiently large conduction heterogeneities, such as large coronary arteries, is activated. In order to understand the mechanism behind LEAP, we have carried out a statistical study of resetting a medium filled by one or more stable spirals (‘rotors’) in a two-dimensional electrophysiological model of cardiac tissue perforated by blood vessels to the resting state (‘defibrillation’). We found the highest success probabilities for this defibrillation for underdrive pacing with periods 10–20 percent larger than the dominant period of the stable rotors in the unperturbed dynamics. If a sufficiently large number pulses is applied and an optimal pacing period chosen, the energy per pulse required for successful defibrillation is about 75–80 percent lower than the energy needed for single-shock defibrillation. Optimal conditions to control and suppress fibrillation based on stable rotors, hence, are similar to the ones in found for the case of an electrophysiological model displaying spatiotemporal chaos (‘electrical turbulence’) in an earlier study (see P Buran et al 2017 Chaos 27 113110). The optimal pacing period is found to increase with increasing strength of the electrical field strength used in the model. The success probability also increases strongly until the fourth or fifth pulse administered, which is strongly correlated to an observed increase of the fraction of re-excitable tissue with each subsequent pulse. Monitoring the fraction of excitable tissue in the model as key quantity of the excitable medium, moreover, enabled us to successfully predict the optimal pacing period for defibrillation.


Introduction
A loss of rhythm and synchronization of the cardiac electrical activity, orchestrating heart contraction and the pumping of blood, is associated with a number of arrhythmias including atrial fibrillation (AF), ventricular tachycardia (VT) and ventricular fibrillation (VF). VF is a particularly dangerous malfunction of the heart, preventing an effective mechanical contraction of the ventricles and causing sudden cardiac death within a few minutes if left untreated. The electrical activity during VF is dominated by multiple stable spirals (rotor) or by multiple wavelets in a spatiotemporally chaotic state [1][2][3][4]. VT on the other hand is often assumed to be generated by one single rotor; VT is not directly lethal, it can, however, degenerate rapidly to VF.
Defibrillation by a strong electrical pulse is the standard way of terminating VF. However, strong shocks are associated with adverse side effects including pain and trauma [5] as well as damage of the myocardium [6]. These adverse effects could be potentially avoided or diminished if VF would be terminated reliably by alternative method using significantly lower energy. Recent experimental studies [7,8] of AF in vitro and in vivo and VF in vitro have shown that a sequence of five or more far-field pulses with stimulation rates close to the arrhythmia cycle length typically require 80%-90% less energy per pulse than needed for defibrillation with a single shock. This method has been named low-energy antifibrillation pacing (LEAP). A similar energy reduction was also found with much faster pacing rates than the cycle length for AF [9,10], while the energy may be even reduced further by an order of magnitude for VT [11,12].
Defibrillation is based on the fact that an electric field depolarizes and hyperpolarizes the tissue close to conductivity heterogeneities [13] which act as virtual electrodes [14]. The strength of this effect depends both on the strength of the electrical field and on the size and shape of the heterogeneities [15][16][17]. Only tissue at major conduction heterogeneities may be activated by low-energy pulses. Therefore, global tissue activation and wave termination often originate only from a few localized activation sites (hot spots). Luther et al [8] suggested that these hot spots are given by large coronary arteries and showed that the size distribution of the heterogeneities follows a power law. The results of Caldwell et al [18], however, did contest such a co-localization of hot spots and major coronary vessels.
The mechanism how LEAP terminate arrhythmias is not well understood yet since experimental methods to visualize the whole three dimensional electrical activity inside the heart are missing. A new promising method which might give deeper insights in the future is the determination of the three dimensional electrical activity over a ultrasound measurement of the mechanical activity of the heart [19]. Current theoretical approaches to understand LEAP often focus on the process of unpinning and removal of a small number of stable spirals [15,[20][21][22]. In these papers, it was demonstrated that spirals that were pinned at larger heterogeneities are best terminated if the LEAP pacing frequency is set to 80%-90% of the rotation frequency of the spiral around the pinning site. This choice allows for an efficient scanning of the phase of the spiral in order to find the vulnerable window for spiral termination [23,24]. Another common interpretation of the mechanism behind LEAP is synchronization, suggested by the observations in experiments [8,25] and simulations [25,26] wherein each pulse gradually entrains during successful LEAP more tissue until the entire tissue is synchronized and no excitation waves are left. In this studies LEAP was successful for overdrive [8,26] as well as for underdrive [8,25] pacing.
In a previous paper [27] we have carried out a statistical study of success probabilities for low energy periodic pacing applied to states of spatiotemporal chaos within exemplary two dimensional electrophysiological models, in which the non-conducting heterogeneities were represented by circles whose size distribution was given by the experimental distribution of radii of coronary arteries reported in reference [8]. Two alternative cell models [28,29] exhibiting similar length scales of electrical turbulence and the correct magnitude of the electric field necessary single pulse defibrillation were studied in detail. In both cases, a near-to-resonant pacing with the dominant period of unperturbed electrical turbulence results in maximal energy reduction already at moderate field strength. For a model featuring a pronounced dominant period, the reduction reached up to more than 98% per pulse compared to single pulse defibrillation. For larger field strengths, moderate underdrive pacing was found to be somewhat more efficient than resonant pacing.
In the present study we use again a simple two dimensional model of cardiac tissue perforated by blood vessels, however, focus on the termination of a regime of stable rotors, using a modified version of the Luo-Rudy model [28] with parameter values from [30]. Therein, the defibrillation success probability of a periodic pacing protocol analogous to the LEAP method was determined starting from a state with a single free spiral, a single pinned spiral or multiple free spirals. The aim of this study is to get a deeper understanding of LEAP for stable spirals and to find the optimal protocol in this case. Therefore we have investigated the defibrillation success rate and other aspects of LEAP. It turned out that the optimal pacing period with respect to the termination probability is typically slightly larger than the dominant period of the stable spirals (=rotation period), i.e. a mild so-called underdrive pacing provides optimal conditions to terminate the activity starting from a state of one or more stable rotating spirals. In addition, a high fraction of excitable tissue right before the application of the pulses was found to be reliable indicator for the success of defibrillation. This knowledge allowed for a computationally much more efficient alternative to estimate the optimal pacing period for LEAP.
In section 2, we introduce the cellular model along with the numerical methods, our choice of initial conditions and for the periodic protocol. Section 3 contains the main results regarding single pulse high-energy defibrillation as a reference and defibrillation by low-energy periodic pacing as the main subject of this study. Section 3 concludes with the demonstration that the success of defibrillation by periodic pacing is strongly correlated with the amount of available excitable tissue, which depends on the time since the previous pulse and the number of pulses applied. We show that this observation allows a reliable prediction of the optimal pacing period with the highest success rate. Section 4 provide a discussion of the results and a short conclusion, respectively. Snapshot for a state with a single spiral, whose cores is located at the center. (c) Snapshot for a state with a single pinned spiral, whose cores is pinned to a circular heterogeneity at the center with a radius of 0.025L. Conductivity heterogeneities are represented by small blue circles. Movies of the simulations are provided in the supplementary material (https://stacks.iop.org/NJP/24/083024/mmedia).

Methods
We have performed extensive numerical simulations for defibrillation protocols with different electrical field strengths, number of pulses and pacing periods for states with multiple free stable spirals, one free stable spiral and one pinned stable spiral. In order to obtain statistically valid results, typically 50 simulation runs with different initial conditions were performed for each configuration. Therefore, we restricted the simulations to the same simple two-dimensional model of homogeneous, isotropic cardiac tissue, used in our previous work [27], in which the heterogeneities are represented by circles whose size distribution is given by the distribution of radii of coronary arteries measured in reference [8] and which bases on the monodomain model [31] in conjunction with a cellular model.

Numerical methods and model parameters
Simulations were carried out on a quadratic two-dimensional equidistant finite-difference grid with no-flux boundary conditions and 1000 × 1000 nodes, resulting in a grid spacing of Δx = 10 −3 L, where L is the system size (side length of the quadratic tissue domain). Conductivity heterogeneities were included into this grid as non-conducting patches by modified Neumann-boundary-conditions [15,17] and the time integrations were performed with a time step of 0.1 ms, which was sufficient to guarantee numerical stability.
In order to keep the results comparable with the ones in our previous work [27], we used the exact same configuration of the simulation domain regarding system size L and distribution of heterogeneities. The system size of the quadratic tissue domain was L = 20 cm, which has then the extension of the order of magnitude of the surface of a (large) mammalian heart. The distribution of the heterogeneities followed the experimental size distribution p(R) ∝ R −2.75 of large coronary arteries given in Luther et al [8]. As the cutoffs of the power law and as the density of all heterogeneities R min = 3 × 10 −4 L, R max = 4 × 10 −3 L and ρ 0 = 1.6 × 10 4 /L 2 were chosen. Note, that the area covered by defects amounts to around 3 percent of the total area used in the simulations in this paper. Heterogeneities with a diameter smaller than the grid spacing (2R < Δx) where set to the size of one pixel. We have checked that the usage of the rather low spatial resolution Δx = 0.001L and also the modeling approach of setting heterogeneities with a diameter smaller than the spatial resolution to the size of one pixel, does not significantly impact the simulation results, see appendix A. We have further checked that different realization of positions and radii resulted in almost identical defibrillation successes. Therefore, the same realization of positions and radii was used in simulations while the initial states were varied unless stated otherwise.
As the cellular electrophysiological model, we used the Luo-Rudy model [28] with parameter values from [30], which exhibit stable spirals. This allowed us to study initial conditions with multiple spirals (associated with fibrillation) and single spirals either free or pinned (which are often associated with tachycardia), see figure 1. The diffusion constant with a value of D = 5.625 × 10 −3 L 2 s was chosen such, that it resulted in a single defibrillation threshold that is consistent with experiments [8,32], see figure 3. Scaled spectra of transmembrane potential V for multiple spirals (blue), a single spiral (red) and a single pinned spiral (green). Simulations were performed for 10 s and with 10 different initial conditions. The displayed graphs are averages of spectra of the transmembrane potential V at all spatial grid points and from all the 10 simulations runs. The spectra are in all threes cases nearly identical and characterized by a narrow peak at f 0 = 10 Hz and further local maxima at multiples of this frequency.

Figure 3.
Single pulse termination probabilities P S for multiple spirals (blue), a single spiral (red) and a single pinned spiral (green) as functions of the applied field strength E and for biphasic (circles) and monophasic (squares) pulse shapes with 10 ms duration. The dotted horizontal line indicates a termination probability of 90% and defines E 90 , i.e., the minimum electric field strength to obtain a termination probability of at least 90%.

Determination of fraction of excitable tissue
The fraction of excitable tissue of a state was calculated as the number of excitable points divided by the total number of points. Excitable points were characterized by having a transmembrane potential lower than the activation threshold (V thresh = −40 mV) and releasing an action potential if the transmembrane potential exceeds this activation threshold. Therefore, additional time integration had been performed, where the transmembrane potential was set to the value of the activation threshold, and it was checked whether the transmembrane potential would achieve a value greater than 0.0 mV.

Defibrillation protocols
State-of-the-art defibrillators execute a biphasic, asymmetric protocol [33]. Nevertheless, we initially tested both monophasic and biphasic defibrillation protocols. For monophasic pulses, a rectangular waveform with a amplitude E and duration of 10 ms was chosen. The waveform for biphasic pulses consisted of two subsequent rectangular parts with identical amplitudes E and opposite field direction. We fixed the total pulse duration to 10 ms and found that 7 ms (forward) and 3 ms (backward) lead to the best defibrillation results, similar to our results found in electrophysiological models exhibiting spatiotemporal chaos [27]. The used protocols contained n identical biphasic pulses with a constant interval and the pacing period T is defined as the time difference between the onset of subsequent pulses.

Generation of initial conditions
A single free spiral wave was initiated using the standard cross-field stimulation technique [34]. We run the simulation for 5 s to allow the spiral(s) to stabilize. Subsequently, we let the simulation run for one more complete rotation of the spiral to obtain the single spiral initial states at different equidistant phases. For the pinned spiral states, we placed one circular heterogeneity with a radius of 0.025L at the center of a single spiral state. Then, the simulation was carried out for 5 more s to allow the spiral to pin and to stabilize. Subsequently, we ran again the simulation for one more complete rotation of the pinned spiral to obtain the . Sequence of snapshots for successful pacing with a period of T/τ 0 = 1.15 and field strength E = 3.5 V cm −1 of a single pinned spiral, whose core is pinned to a circular heterogeneity at the center. The snapshots show the transmembrane potential right before, during, and 80 ms after the first, second, third, fourth and fifth pulse. (a) State right before a pulse. (b) The pulses induce excitation fronts at the big heterogeneities, which excite quickly nearly the entire excitable tissue, preventing the spirals to spread. (c) As a result of this, big spirals fall apart into smaller spirals while small spirals might get completely eliminated. Overall the by each pulse excited tissue fraction increases with every pulse until all excitation fronts are terminated. pinned spiral initial states at different equidistant phases. In order to get the multiple spirals initial states, we first cut from five randomly chosen single spiral states a quadratic section. Each section had a side length of 500 nodes (half of the system size) and its center was located in the middle of the single spiral, randomly additionally shifted up to 25 nodes in each direction. Then we placed in every quarter of a quiescent state one and at a random position the fifth of these five sections of cut spirals. Finally we run this generated states for 1 s to stabilize. In all cases, we generated 50 different initial states.

Performed simulations
In order to obtain statistically valid results, typically N = 50 simulation runs with different initial conditions were performed for each parameter configuration and type of stable-spiral-state (multiple free spirals, a single free spiral and a single pinned spiral). Already for the systematic analysis of periodic pacing dynamics of stable spirals in section 3.4 alone, 78 000 simulations runs with a sequence of 5 pulses were performed for: For the entire study a total of about 100 000 simulation runs were performed. Figure 2 shows the spectra for states with multiple spirals, a single spiral and a single pinned spiral. All three spectra are nearly identical and do not seem to depend on the number of spirals nor on the fact whether a spiral is pinned or not. They are all three characterized by a very narrow peak at f 0 = 10 Hz or τ 0 = 100 ms respectively corresponding to the rotation period of the spirals and further local maxima are at multiples of this frequency. The characteristic length scale of the unperturbed patterns is given by the wavelength of the spiral(s) that amounts to around 5.5 cm. This scale is almost two orders of magnitude bigger than the radius of the largest defect which equals 0.08 cm. As a result the defects in the medium do not change propagation properties substantially.

Single pulse defibrillation
We determined the single pulse termination probability P S , which is defined as the probability that a fibrillation state is converted into a quiescent state by a single defibrillation pulse, as a function of field strength E. We computed P S as the fraction of successful termination events by running N = 50 simulations (see supplementary material) with different initial states for each parameter configuration, see figure 3. The statistic error of P S is given by We have checked both monophasic and biphasic protocols, determining that biphasic (3 ms forward, 7 ms backward) are more effective than monophasic (10 ms forward) pulses, see figure 3. This is in line with earlier findings [35] and our previous results for different models exhibiting spatio-temporal chaos (=electrical turbulence) [27]. The probability to terminate a single free or pinned spiral is found to be almost identical than the one for terminating multiple spirals. Furthermore, the minimum electric field strength to reach a termination probability of at least 90% is nearly identical independent of the initial state. We found E mono 90 = 9.5 V cm −1 for monophasic and E bi 90 = 5.9 V cm −1 for biphasic pulses. Our results roughly compare with systematic experiments with monophasic pulses applied to over 200 fibrillation episodes of canine hearts in vivo. While defibrillation is hopeless for field strengths below 2.9 V cm −1 , the probability to defibrillate increases for higher values and is guaranteed for field strength above 8.5 V cm −1 [32].

Periodic anti-fibrillation pacing
We employed a sequence of biphasic pulses to mimic the experimental LEAP approach, since they require a lower field strength for defibrillation compared to single monophasic pulses, see section 3.2.
Taking a first look on successful LEAP of a single pinned spiral (figure 4, we observe on the one hand that the spiral splits first up into smaller spirals until they get consecutively eliminated and on the other hand that the by each pulse excited tissue fraction increases with every pulse. Furthermore, we could not observe an unpinning process of the pinned spiral nor a shifting of the free spirals to the boundary. Similar observations could be made for multiple or a single free spiral, see provided videos in the supplementary material. This indicates already that neither the number of spirals nor the fact whether a spiral is pinned or not plays an essential role in the success of LEAP.
The choice of the pacing period T turns out to be crucial for LEAP success, see figure 5. The maximal termination probability P L is found for underdrive pacing with a period of T max = 115 ms, which is 15% higher than the dominant period τ 0 = 100 ms corresponding to the rotation period of stable spirals in the medium, see figure 2. Further local maxima are found for T max /2 and 2T max albeit with considerable lower termination probabilities. The termination probability P L increased substantially with the number of pulses n. It is worth noting, that this picture is remarkable close to the success probabilities for defibrillation found in our earlier study using models exhibiting electrical turbulence [27], though therein shape of the observed peaks was considerably more smeared out.
Choosing the optimal pacing period T max , the minimal field strengths E L 90 for a defibrillation success probability of at least 90% increases monotonously with the number of pulses n, see the colored dotted lines in figure 6. We found for multiple spirals values of E L 90 ≈ 3.5 V cm −1 for five, E L 90 ≈ 2.8 V cm −1 for ten and E L 90 ≈ 2.7 V cm −1 for fifteen pulses. Further substantial reductions for more pulses are unlikely as the difference between fifteen and ten pulses is already very small. The overall energy delivered by a periodic pacing protocol is proportional to the number of pulses n and the square of the electric field strength E. Compared to the field strength E bi 90 = 5.9 V cm −1 required for a single defibrillation pulse, the total required Energy for LEAP is always higher then for a single pulse. However there is a reduction of dissipated energy per pulse, which is for multiple spirals about 65% for five, 77% for ten and 79% for fifteen pulses.
The required energies for a single free or pinned stable spiral are for more then three pulses almost identical and only a little bit higher then the one required for multiple spirals. Hence, neither the number of spirals nor the fact whether a spiral is pinned or not seems to play an essential role for the success of periodic pacing.

Analysis of periodic pacing dynamics
In the following we will only look at periodic pacing protocols with up to n = 5 pulses as the qualitative behavior does not change much for higher number of pulses. If not mentioned otherwise, averages and the termination probability P L are computed by running N = 50 simulations with different initial states for each parameter configuration. Moreover we will only show the results for multiple spirals as the results for a single free or pinned spiral were similar.

Optimal pacing period
The range of successful pacing periods T, for which defibrillation is found in the simulations, becomes wider for high electric field strengths E, see figure 7(a). We define the lower bound T − L and upper bound T + L of this range as the two positions with 80% of the maximal possible termination probability: This allows us to study pacing protocols in regimes with low termination probabilities as well, compared to take the two positions at a fixed termination probability. Taking the global maximum T max is not a good characterization for the optimal pacing period since it might not be located at the center of the peak. Therefore we define the optimal period T opt as the mean of the two bounds T − L and T + L : The range of pacing periods T with more than 80% success rate starts at a minimal field strength, where the lower and upper bounds intersect, and widens monotonously with increasing field strength E, see figure 7(b). For low E, the lower and upper bounds both increase at roughly equal rate. However, for high E the lower bound remains almost constant, while the upper bound continuous to increase. As a result of that, the optimal period remains first constant and increases gradually for high E. For increasing number of pulses n the region of successful pacing periods becomes wider, while the minimal field strength decreases.

Experimental setup for an efficient estimation of the optimal pacing period
One can observe that during successful periodic pacing big spirals disintegrate into smaller spirals while small spirals get successively eliminated with each pulse until everything is terminated, see figure 4. Hence, a progressively decrease of the total size of the spirals, given by the complement of the fraction of excitable tissue F Exc seems to be the crucial factor for the success of periodic pacing. Successful pacing protocols are characterized by high fractions of excitable tissue F Exc right before the pulses. Achieving a maximal fraction of excitable tissue F Exc right before a pulse is a trade-off between pacing fast enough in order to prevent already recovered tissue to get excited by remaining excitation fronts and pacing slow enough, giving the tissue enough time to fully recover. As shown in figure 8(a), the range of pacing periods T with a high average fraction of excitable tissue right before a pulse F Exc becomes for an increasing pulse number n and electric field strength E wider towards higher values of T. That is, since the total size of the remaining excitation waves decreases during successful pacing, see figure 4(a), slowing down the excitation of excitable tissue. This widening of the range of pacing periods T with a high average fraction of excitable tissue right before a pulse F Exc towards higher values of T looks similar to the one for the range of high termination probabilities P L , see figure 7(a). As shown in figure 8(b), the average fraction of excitable tissue right before a pulse F Exc is highly correlated with the termination probability P L of a pulse, at least for pulse numbers n 3. A good estimation of the optimal pacing period T opt , which is optimal in the sense of a high termination probability P L , should be therefore given by a optimal pacing Figure 9. (a) Optimal pacing period T opt (black) with respect to a high termination probability P L and optimal pacing period T Exc opt (red) with respect to a high average of the by pulse excited tissue fraction F Exc as a function of field strength E and different pulse numbers n. T Exc opt is a good estimation of T opt for high E. It lies for n = 3 and n = 4 and field strength lower then 4 V cm −1 and 3.4 V cm −1 respectively outside the regime of successful pacing periods, however, coincides for n 5 almost completely with T opt . (b) Same as in (a), however, F Exc is here determined by only one randomly chosen initial condition and not as an average for every parameter configuration (E, T, n). Nevertheless, the impact on T Exc opt are not significant and T Exc opt keeps a nearly equally good estimation of T opt compared to the one from (a) for an averaged F Exc . period with respect to high F Exc . We define this pacing period T Exc opt analogous to T opt as the mean of the lower and upper bound T ± Exc of pacing periods T achieving a F Exc of at least 80% of the maximal possible F Exc : As shown in figure 9(a), T Exc opt is a good estimation of the optimal pacing period T opt and coincides for pulse numbers n 5 almost completely with T opt . This also applies for lower n in the regime of high field strength E and T Exc opt lies for this field strength E, where the correlation coefficient of P L and F Exc is greater then 0.8, see figure 8(c), still in the regime of successful pacing periods.
Measuring the fraction of excitable tissue right before a pulse F Exc is harder then to determine whether a pulse is successful or not, especially in a real experiment. Therefore, it would not be worth to determine the optimal period T opt by the use of F Exc instead of the success probability P L . However the fraction of excitable tissue right before a pulse F Exc of a single pulse can be measured directly, while P L can be only determined as the mean of successful (0) and unsuccessful (1) events. This might allow us to determine an accurate value for F Exc with less simulations or experiments and determining F Exc by only one randomly chosen initial condition for every parameter configuration has already no significant impact on T Exc opt , as shown in figure 9(b). T Exc opt still coincides for pulse numbers n 5 quiet good with T opt and is only a little bit noisy.

Conclusion and discussion
We have carried out a statistical study of simulations of periodic pacing for an electrophysiological model exhibiting stable spirals within a two-dimensional model of cardiac tissue perforated by blood vessels. It turns out that underdrive pacing (=pacing with a period lower than the spiral rotation period) results in maximal energy reduction and that neither the number of spirals nor the fact whether a spiral is pinned or not plays an essential role in the defibrillation success of periodic pacing. A protocol with ten pulses and a pacing period 15% higher than the dominant period requires 77% less energy per pulse than a single biphasic shock for successful defibrillation. We observed that successful pacing protocols are characterized by repeated pulsing at times when the fraction of excitable tissue is high. This fraction of excitable tissue right before a pulse is found to be highly correlated with the termination probability allowing an efficient alternative to estimate the optimal pacing period to achieve defibrillation. The appropriate choice of length scales is an important issue. We have chosen a scaling based on the experimental observation regarding the sizes of the heterogeneities in [8] and picked a value of the conductivity such that the single pulse defibrillation threshold reproduces experimentally found values [8,32]. Furthermore, the size of the system was picked sufficiently large to generate initial states not only with a single free or pinned spiral but also with around ten coexisting spirals. While the initial dynamics after the start of pacing depended strongly on the initial conditions the success probability after a sufficiently large number of pulses was similar for all chosen initial conditions. We found that defibrillation by periodic pacing works most efficiently with underdrive pacing for stable spirals compared to our results for spatiotemporal chaos [27] where resonant pacing was favorable at low field strength and underdrive pacing was only superior for stronger electrical field. An increase of the optimal pacing period with the field strength is also observable for stable spirals. This is caused by an expansion and not by a shift of the regime of successful pacing towards higher values of the period T and is in line to our earlier observation for defibrillation by periodic pacing in models exhibits spatiotemporal chaotic electric activity [27]. A protocol with 10 pulses and a pacing period 15% higher than the dominant period requires 77% less energy per pulse than a single biphasic shock for successful defibrillation, which corresponds with the energy reductions between 75% and 90%, found for successful underdrive pacing of canine hearts [8].
Current theoretical approaches to understand LEAP for stable spirals usually focus on the process of unpinning and the removal of free stable spirals by shifting the spiral tips to the boundary of the heart [15,16,20,21]. These theories predict also optimal success for underdrive pacing with periods between 110% and 120%, [23,24]. However, our observation, that neither the initial number of spirals nor the fact whether a spiral is pinned or not plays an essential role in the success of periodic pacing, see figures 5 and 6 point towards different essential ingredients. Furthermore we found in our simulations that the process of unpinning and shifting of the spirals to the boundary of the heart did not play an essential role, see figure 4 and provided videos in the supplementary material. In contrast, we observed that during successful pacing big spirals disintegrate into smaller spirals while small spirals get successively eliminated with each pulse until all activity is terminated, see figure 4.
Therefore, a progressive decrease of the total size of the spirals in parallel to an increase of the fraction of excitable tissue F Exc turned out to be the crucial factor for the success of periodic pacing. In other words, successful pacing protocols are characterized by pulsing at times when the fraction of excitable tissue F Exc is high. This picture is consistent with observations, made in experimental [8,25] and numerical studies [25,26], which interpret repeat pacing as a synchronization process of the phases of the electrical activity in the tissue. The subsequent increase in the total fraction of excitable tissue with each pulse in the periodic pacing leads to an increasing synchronisation of local excitation cycles, which is, however, depending on the simultaneous trigger provided by the pacing and therefore different from the more frequently studied mutual synchronisation occurring in oscillatory media.
The observation that the termination probability P L of a pulse is highly correlated with the fraction of excitable tissue right before a pulse F Exc for pulse numbers n 3, see figure 8(b) allows us to determine a criterion for the optimal pacing period T opt (optimality here refers a high termination probability P L ). This done by finding the pacing period that maximizes F Exc , see figure 9(a). The big advantage of this estimation is the fact, that the excitable fraction F Exc is easily accessible and could be measured directly in experiments, whereas the termination probability P L can be only determined by counting successful (0) and unsuccessful (1) events requiring therefore a large amount of independent measurements. In our context, determining F Exc by only one simulation for every parameter configuration provides a very good approximation of the optimal pacing period, see figure 9(b).
Finally, it is worth pointing out that the results in this study of an electrophysiological model that represents an excitable medium wherein one or more stable spirals can coexist in the absence of pacing is in many ways very similar to the results that we found earlier for two electrophysiological models exhibiting spatiotemporal chaos (=electrical turbulence) [27]. Therefore, the mechanistic considerations described in the last part of the results section above may generalize to these qualitatively different dynamics.
As described in section 2.1, the heterogeneities in our composite model of two-dimensional homogeneous isotropic tissue represent non-conducting blood vessels that have circular shape. The size distribution of these circles is assumed to follow a power law, that was reported by Luther et al [8] for the vasculature of the heart muscle: The density of activated heterogeneities for a given distribution of sizes (respectively radii) is obtained as: In all our simulations in the results section above, the smallest considered radius was R min = 3 × 10 −4 L, the largest considered radius was R max = 4 × 10 −3 L, the exponent was α = −2.75 and the density of all heterogeneities was ρ 0 = 1.6 × 10 4 , where L is the system length. The heterogeneity radii drawn from this distribution are close to the spatial resolution Δx = 0.001L, which was used for most of the simulations in this study. Heterogeneities with a diameter even smaller than the spatial resolution (2R < Δx) were set to the size of one pixel. Such limitations in resolving circular heterogeneities are expected to mainly affect the activation behavior of the heterogeneities, since it may depend strongly on the size and shape of a heterogeneity used in the simulations [15][16][17].
Resolution-related artifacts are clearly visible for the dependency of the minimal radius R min (E) that is necessary for a heterogeneity to get activated by a pulse on the applied electric field strength in figure A1. The minimal activation radius R min (E) decreases step-wise and not smoothly for the standard spatial resolution Δx = 0.001L (green), staying constant within a wide range of field strength E. These resolution-related artifacts do not occur anymore for R min (E), when the spatial resolution Δx = 0.0001L (yellow) is chosen to be ten times higher.
Resolution-related artifacts may hence be expected to appear also for the density of activated heterogeneities ρ according to relation (A.2). However, the relation (A.2) does not take into account further resolution-independent factors like the current excitability of the surrounding tissue of a heterogeneity or whether there are other heterogeneities nearby that support or counteract the depolarization of the surrounding tissue. Such factors have also a decisive influence on whether a heterogeneity gets activated, and thus can compensate or at least mitigate the resolution-related factors considered in (A.2). Indeed, for the activation density ρ determined directly as the average quotient of the number of activated heterogeneities after a pulse and the area size of excitable tissue right before the pulse, resolution-related artifacts are almost not observable, see figure A2. Moreover, figure A2 shows that the activation density ρ is independent from the fact whether it is determined in simulations where all heterogeneities with a diameter smaller than the grid spacing (R < Δx/2 = 0.0005L) were set to the size of one pixel (blue), or whether it was determined in simulations where only heterogeneities with R Δx/2 = 0.0005L were considered (cyan). In other words, the slight distortion of the size distribution of heterogeneties that stems from the relatively coarse discretization used in our simulation described above does not affect the activation density ρ. Figure A1. Minimal radius R min (E) that is necessary for a heterogeneity to get activated by a pulse as a function of the electric field strength E of the pulse. The radii R min (E) were determined in single heterogeneity simulations, where it was checked for each field strength E and a wide range of radii whether a circular heterogeneity that is placed in the center of the simulation domain gets activated or not by a pulse. For the spatial resolution Δx = 0.001L (green), which is used by most of the simulations in this study, resolution-related artifacts are clearly visible: R min (E) decreases step-wise and not smoothly, staying constant within a wide range of field strength E, in contrast to R min (E) determined for a ten times higher spatial resolution of 0.0001L (yellow), where those resolution-related artifacts do not occur anymore. Figure A2. Density of activated heterogeneities ρ, as a function of the electric field strength E. ρ determined directly in simulations of each initial state as the average quotient of the number of activated heterogeneities after a pulse and the area size of excitable tissue right before the pulse, is plotted in blue for simulations where all heterogeneities with a diameter smaller than the grid spacing (R < Δx/2 = 0.0005L) were set to the size of one pixel and in cyan for simulations where only heterogeneities with R Δx/2 = 0.0005L were considered and placed on the simulation domain. Both curves are almost equal for field strength E below the single pulse defibrillation threshold (vertical dotted line) indicating that the chosen approach of setting heterogeneities with R < 0.0005L to the size of one pixel should not significantly impact the simulation results for LEAP. Furthermore, resolution-related artifacts due to the low spatial resolution of Δx = 0.001L, which occur for the minimal radius that is necessary for a heterogeneity to get activated R min (E) (see figure A1) and which could also be expected to occur for ρ (due to relation (A.2)), are not observable: the course of both directly measured ρ (cyan and blue) is smooth and does not increase step-wise as one would expect at first, estimating ρ with equation (A.2) via the minimal radius R min (E) (green). Instead, the course is more similar to the estimation of ρ for a ten times higher spatial resolution of 0.0001L (yellow), for which resolution-related artifacts of the underlying minimal radius R min (E) do not occur anymore.
Furthermore it can be seen, that the approach how heterogeneities with a diameter smaller than the grid spacing are modeled has no significant impact on the density of activated heterogeneities ρ, at least not for field strength E below the single pulse defibrillation threshold, which were mainly used in this work. The course of ρ determined in simulations where all heterogeneities with a diameter smaller than the grid spacing were set to the size of one pixel (blue) is almost equal to the course of ρ determined in simulations where, in contrast, all these small heterogeneities were removed (cyan).
Overall, it can be therefore assumed that the usage of the spatial resolution Δx = 0.001L should not significantly impact the simulation results in this work, since the density of activated heterogeneities ρ, which mainly characterizes the activation behavior of the heterogeneities, is neither significantly affected by the limitations of resolving heterogeneities correctly, nor by the modeling approach of setting heterogeneities with a diameter smaller than the grid spacing to the size of one pixel.