The role of pulse timing in cardiac defibrillation

Life-threatening cardiac arrhythmias require immediate defibrillation. For state-of-the-art shock treatments, a high field strength is required to achieve a sufficient success rate for terminating the complex spiral wave (rotor) dynamics underlying cardiac fibrillation. However, such high energy shocks have many adverse side effects due to the large electric currents applied. In this study, we show, using 2D simulations based on the Fenton-Karma model, that also pulses of relatively low energy may terminate the chaotic activity if applied at the right moment in time. In our simplified model for defibrillation, complex spiral waves are terminated by local perturbations corresponding to conductance heterogeneities acting as virtual electrodes in the presence of an external electric field. We demonstrate that time series of the success rate for low energy shocks exhibit pronounced peaks which correspond to short intervals in time during which perturbations aiming at terminating the chaotic fibrillation state are (much) more successful. Thus, the low energy shock regime, although yielding very low temporal average success rates, exhibits moments in time for which success rates are significantly higher than the average value shown in dose-response curves. This feature might be exploited in future defibrillation protocols for achieving high termination success rates with low or medium pulse energies.

1 Introduction

Sudden cardiac death, in many cases due to ventricular fibrillation (VF) Koplan and Stevenson (2009), causes around 700,000 deaths per year in Europe alone, and hence, remains one of the major public health issues Chugh (2017). In practice, such cardiac arrhythmias are most successfully terminated by applying external electric shocks to the heart, referred to as defibrillation. Demonstrated both experimentally Davy et al. (1987); Kwaku and Dillon (1996) and numerically Boegli et al. (2000); Bragard et al. (2013); Lilienkamp and Parlitz (2020), the dependence of the defibrillation success rate on the strength of the shock applied, which is called dose-response curve (DRC), is given by a sigmoid function. However, high energy shocks, being promising in terms of successful defibrillation, cause severe pain and tissue damage to the patient Holden et al. (2019). Another adverse side effect is that postshock arrythmias are much more probable for high defibrillation field strengths Cates et al. (1994); Fast and Cheek (2002). A major challenge in medical practice is thus, to minimize these tremendous side effects by significantly reducing the shock energy while keeping the probability for successful defibrillation high. Therefore, different types of approaches for energy reduction have been proposed, including a sequential application of low energy pulses Luther et al. (2011) or multi-stage electrotherapy Janardhan et al. (2014). In contrast, the current study points the way to a possible improvement in low-energy single-shock therapy. Using computational methods, we show that the success rate of external shocks may strongly fluctuate, a feature that might be exploited for the defibrillation of arrhythmic dynamics of the heart using low energy shocks. In our study, we used the Fenton-Karma model Fenton and Karma (1998), describing the electrical excitation dynamics of cardiac cells in a 2D computational domain and simulated chaotic spiral waves of cardiac excitation corresponding to fibrillation in a real heart. By means of the concept of virtual electrodes Trayanova et al. (1998); Pumir and Krinsky (1999); Bittihn et al. (2008), we then perturbed the spiral wave dynamics, aiming at terminating it using pulses of low energy. To conduct this study, we generated several trajectories of chaotic spiral wave dynamics and studied the time series of the success rate for terminating the complex dynamics by applying such perturbations. In this way, optimal times for the termination of cardiac arrhythmias were found. In contrast, the rather small success rate for low-energy shocks observed in the DRCs is a consequence of applying shocks at random times and averaging over these attempts. In terms of real-life defibrillation, this would mean that, if external low-energy shocks are applied at the right moment in time, defibrillation success rates might be significantly higher than when applying them at random times. By preventing the aforementioned negative side effects of high-energy defibrillation strengths, this could lead to a major improvement in medical application. The paper is structured as follows. In Section 2, we introduce the Fenton-Karma model, its implementation and how defibrillation is simulated. Section 3 reproduces the DRC qualitatively, which is then extended to a study of how the values of the success rate for fixed perturbation parameters are distributed temporally, i.e., within the time series. Afterwards, we perform a robustness study with respect to the fixed perturbation control parameters chosen herein and show how the waiting times for high success rates are distributed. Finally, in Section 4, we summarize the results obtained and critically discuss their generality and applicability to real-life defibrillation.

2 Methods

In this section, we will describe the model used for the study, as well as how it is implemented numerically. We will furthermore show how the complexity of chaotic cardiac dynamics (corresponding to fibrillation) is quantified in terms of phase singularities and explain how defibrillation is realized in our simulations.

2.1 The Fenton-Karma model and its numerical implementation

For the study of spatiotemporally chaotic spiral waves of cardiac excitation representing fibrillation, we choose a 2D implementation of the Fenton-Karma model Fenton and Karma (1998). Its dynamics is governed by the following reaction-diffusion equations describing the temporal evolution of the membrane potential V_m and the gating variables, v and w:

$\frac{\partial V_{\text{m}}}{\partial t}=\nabla \cdot \tilde{\mathit{D}}\nabla V_{\text{m}}-\left[I_{\text{fi}}\left(V_{\text{m}},v\right)+I_{\text{so}}\left(V_{\text{m}}\right)+I_{\text{si}}\left(V_{\text{m}},w\right)\right]/C_{\text{m}},(1)$

$\frac{\partial v}{\partial t}=\mathrm{\Theta }\left(u_{\text{c}}-V_{\text{m}}\right)\left(1-v\right)\left[\frac{\mathrm{\Theta }\left(V_{\text{m}}-u_{\mathit{v}}\right)}{{\tau }_{\mathit{v}1}^{-}}+\frac{\mathrm{\Theta }\left({\mathit{u}}_{\mathit{v}}-{\mathit{V}}_{\text{m}}\right)}{{\tau }_{\mathit{v}2}^{-}}\right]-\mathrm{\Theta }\left(V_{\text{m}}-u_{\text{c}}\right)\frac{v}{{\tau }_v^{+}},(2)$

$\frac{\partial w}{\partial t}=\mathrm{\Theta }\left(u_{\text{c}}-V_{\text{m}}\right)\frac{1-w}{{\tau }_w^{-}}-\mathrm{\Theta }\left(V_{\text{m}}-u_{\text{c}}\right)\frac{w}{{\tau }_w^{+}},(3)$

where $\tilde{\mathit{D}}$ denotes the diffusion tensor. For our work, we replaced the standard Heaviside functions Θ(x) by a continuous approximation in order to ensure differentiability Bittihn (2015):

$\mathrm{\Theta }\left(\mathit{x}\right)\to \frac{1}{2}\left(1+\tanh \left[k_{2}x\right]\right),(4)$

with the smoothing parameter chosen to be k₂ = 10. The fast inward (fi), slow outward (so) and slow inward (si) ion currents appearing in Eq. 1 read:

$I_{\text{fi}}\left(V_{\text{m}},v\right)=-\frac{v}{{\mathrm{\tau }}_{\mathrm{d}}}\mathrm{\Theta }\left(V_{\text{m}}-u_{\text{c}}\right)\left(1-V_{\text{m}}\right)\left(V_{\text{m}}-u_{\text{c}}\right),(5)$

$I_{\text{so}}\left(V_{\text{m}}\right)=\frac{V_{\text{m}}}{{\tau }_{\text{o}}}\mathrm{\Theta }\left(u_{\text{c}}-V_{\text{m}}\right)+\frac{1}{{\tau }_r}\mathrm{\Theta }\left(V_{\text{m}}-u_{\text{c}}\right),(6)$

$I_{\text{si}}\left(V_{\text{m}},w\right)=-\frac{w}{2{\tau }_{\text{si}}}\left(1+\tanh \left[k\left(V_{\text{m}}-u_{\text{c}}^{\text{si}}\right)\right]\right).(7)$

Note that the Heaviside functions Θ occurring here are again approximated by the hyperbolic tangent function (Eq. 4). The parameter set used is given in Table 1.

TABLE 1

TABLE 1. Fenton-Karma model parameter set used here Fenton et al. (2002).

We chose a 2-dimensional computational domain, $\mathcal{D}$, with N_x × N_y = 100 × 100 grid points and an equidistant spacing of h_x = h_y = 1 mm. Furthermore, we assumed a constant, homogeneous and isotropic diffusion given by $\tilde{\mathit{D}}=D=0.2{\mathrm{m}\mathrm{m}}^2/\mathrm{m}\mathrm{s}$, which reduces the first term of Eq. 1 to the Laplace operator Δ acting on the membrane potential V_m. The Laplace operator was discretized via a finite differences nine-point stencil Bittihn (2015), while the time-wise discretization was realized with the explicit Euler algorithm Press et al. (2007), using an integration time stepping of dt = 0.2 ms. As boundary condition for the computational domain $\mathcal{D}$, we chose no-flux boundary conditions Bittihn (2015) which were realized by extending $\partial \mathcal{D}$ with ghost points. For generating the spiral wave dynamics investigated, we used a modification of the cross pacing protocol, as described in the Supplementary Material. In order to have a measure for whether a state is governed by (chaotic) spiral waves, we evaluated the phase singularities on $\mathcal{D}$ Iyer and Gray (2001):

$n^{\text{top}}={\frac{1}{2\pi }\mathit{\oint }}_{\mathcal{C}}\nabla \theta \text{d}\mathit{l},(8)$

where $\mathcal{C}$ is a closed integration path along eight grid points around each point of reference and θ is the phase characterizing the dynamics at each grid point Datseris and Parlitz (2022). A more detailed description of the phase, its calculation as well as the parametrization of the integration path is given in the Supplementary Material. Phase singularities are associated with the organizing centers of the spiral wave dynamics of cardiac excitation. Equation 8 yields n^top = 0 if there is no phase singularity enclosed by $\mathcal{C}$ and n^top = ± 1 if there is one enclosed. If at least one phase singularity is detected, we consider the dynamics to be dominated by spiral waves of excitation and thus, to represent arrhythmic dynamics of the heart.

We now want to briefly introduce the concept of virtual electrodes which will be used to simulate the defibrillation process. The heart exhibits heterogeneities of different sizes and configurations, caused, e.g., by blood vessels, discontinuities between sheets of fibers and bundles, or collagen existent in the extracellular space, resulting in a change in tissue conductivity Trayanova et al. (1998). When applying an external electrical field for a sufficiently long time, ions in the tissue are exposed to a Lorentz force which, in turn, shifts them such that they assemble at regions of differing conductivities. If the resulting hyper- and depolarization exceeds a certain threshold, an action potential is triggered at the heterogeneities, resulting in excitation waves Pumir and Krinsky (1999); Bittihn et al. (2008); Lilienkamp (2018). In this sense, in the presence of an electric field, the inhomogeneities act as virtual electrodes. The activation threshold of virtual electrodes depends on the shape (curvature) and size Bittihn et al. (2008); Bezekci et al. (2015) of the heterogeneities and the higher the external field strength the more heterogeneities emit a local excitation wave. This relationship between field strength and the number of activated virtual electrodes is exploited in the current study, where the number of local perturbation sites N_pert is used as a control parameter representing the applied shock energy.

The perturbation sites (representing heterogeneities) were chosen with a fixed size of 2 × 2 grid points and randomly distributed in space without overlap. The effect of the external electric field is approximated by an instantaneous increase of the membrane potential at the perturbation sites

${\left(V_{\text{m},i}^{\text{pert}}\right)}^{2\times 2}={\left(V_{\text{m},i}^{\text{curr}}\right)}^{2\times 2}+{\left(\mathrm{\Delta }{V}_{\text{m},i}\right)}^{2\times 2}\forall i\in \left\{1,\dots ,N_{\text{pert}}\right\}(9)$

with an amplitude of ${(\mathrm{\Delta }{V}_{\text{m}})}^{2\times 2}=0.5$ a.u. The number of perturbation sites N_pert approximately corresponds to the fixed external field strength applied to the cardiac tissue. Finally, as a measure for whether the spiral wave dynamics was successfully terminated, i.e., whether the defibrillation attempt was successful, we checked whether 500 ms after the perturbation (no) phase singularities still existed in the computational domain by evaluating Eq. 8. The absence of phase singularities necessarily leads to convergence of the system towards the steady state, i.e., the electrical excitation dies out on $\mathcal{D}$, such that the regular sinus rhythm could be restored by the sinus node (which is, however, not part of the current model).

3 Results

In what follows, we will summarize our results. First of all, we reproduced the typical sigmoid shape of the dose-response curve for the in silico defibrillation attempts. We then analyzed time series of the termination success rates with a focus on their peak structure, followed by a robustness study of the rather specific choice of perturbation parameters made here.

3.1 Distribution of defibrillation success rates for fixed shock doses

First, we compute the dose-response curve (DRC), describing the dependence of the success rate on the amount of perturbation sites N_pert used to terminate the chaotic excitation dynamics and show that the success rate values are broadly distributed around their mean values. Note that, by making use of the concept of virtual electrodes, N_pert is here considered to directly correlate with the electric field strength applied to the arrhythmic heart. The defibrillation and the determination of whether the attempt was successful was realized as described in Section 2. As a start, we show an example for an unsuccessful perturbation attempt and an example for a successful one (defibrillation). The failed attempt using N_pert = 500 is shown in Figure 1, which also illustrates that the size of the spatial domain is larger than twice the typical diameter of the spiral waves. If, however, not a single phase singularity is detected anymore on $\mathcal{D}$ after 500 ms, the perturbation is considered to be successful like in the example shown in Figure 2 where phase singularities are already absent 30 ms after applying the shock at, in this case, N_pert = 1,000 perturbation sites. In general, a very large number of perturbation sites corresponds to a high shock strength applied to the heart which depolarizes the whole tissue. As a consequence, it is impossible for the spiral wave dynamics to further evolve due to the refractory period not leaving any space for them to propagate.

FIGURE 1

FIGURE 1. An example for an unsuccessful termination attempt for chaotic cardiac excitation using N_pert = 500 randomly distributed perturbation sites (visible as yellow dots immediately after the shock at t = 0 ms). White circles indicate phase singularities. At t = 0 ms the phase singularities were computed immediately before the perturbation. After 500 ms, the state still exhibits phase singularities and since they are considered to be the organizing centers of the spiral wave dynamics, they clearly indicate that the termination attempt failed.

FIGURE 2

FIGURE 2. Successful termination of spiral wave dynamics of cardiac excitation by applying N_pert = 1,000 local stimuli to the same state as shown in Figure 1. White circles indicate phase singularities and were computed for t = 0 ms before the perturbation was applied. Even though the number of phase singularities increases shortly after applying the perturbations (not shown here), they died out already 30 ms after applying the perturbations and the remaining excitation converges quickly to the steady state as shown exemplary at t = 100 ms.

For the DRC, we used 30 different fixed numbers of perturbation sites, N_pert ∈ {50, 100, 150, … , 1500}. For each of these fixed numbers, we randomly generated 50 different spatial configurations and applied them independently to 20 different chaotic states (CS). Thus, for the DRC, we have a total amount of

$\left(\text{no.\hspace{0.17em}of\hspace{0.17em}CS}\right)\times \left(\text{no.\hspace{0.17em}of\hspace{0.17em}}{N}_{\text{pert}}\right)\times \left(\text{no.\hspace{0.17em}of\hspace{0.17em}conf.}\right)=20\times 30\times 50=3\times 10^4$

termination attempts. The corresponding relation between the success rate and the number of perturbation sites applied, S(N_pert), is shown in Figure 3. The shape is indeed sigmoidal and thus qualitatively reproduces the DRCs from experimental Davy et al. (1987); Kwaku and Dillon (1996) and numerical Boegli et al. (2000); Bragard et al. (2013); Lilienkamp and Parlitz (2020) studies.

FIGURE 3

FIGURE 3. Dose-response curve. The dependence of the probability S of a successful termination of the chaotic dynamics on the number of perturbation sites, N_pert, corresponding to the strength of the applied shocks. Each point corresponds to 20 × 50 = 1,000 simulations. The simoid fit S = 100/(1 + c exp (−kN_pert)) is obtained using nonlinear least squares (NLS) to determine the parameters c = 5,632 and k = 0.01106. The inflection point, as well as the average success rate of 4.3 % for the case of using N_pert = 500 perturbation sites are indicated, the latter being referred to in the rest of this study.

Since, due to the averaging over each termination attempt, the representation of the DRC in Figure 3 gives only limited insight into the actual values for the success rate occurring for each fixed N_pert, we also investigated the distribution of S(N_pert) for all fixed N_pert for the same spacing. The resulting success rate distribution DRC is given in Figure 4. It is composed of ten times more defibrillation attempts compared to the averaged DRC shown in Figure 3:

$\left(\text{no.\hspace{0.17em}of\hspace{0.17em}CS}\right)\times \left(\text{no.\hspace{0.17em}of\hspace{0.17em}}{N}_{\text{pert}}\right)\times \left(\text{no.\hspace{0.17em}of\hspace{0.17em}conf.}\right)=100\times 30\times 100=3\times 10^5.$

FIGURE 4

FIGURE 4. Distributions of success rates (color coded) for different numbers of activated perturbation sites, N_pert. Note that the colorbar, indicating the occurrence of a success rate range, is chosen logarithmically and normalized to the most frequently occurring success rate range per N_pert.

As one can see, perturbations with very low or very high N_pert result in narrow distributions of success rate values. Within intermediate ranges of N_pert, on the other hand, especially around the inflection point of the DRC, we have a large range of success rates S(N_pert). This, in turn, implies that the success of a termination attempt is not only dependent on the number of perturbation sites, N_pert, but also on the state (point in time) to which the stimuli are applied.

3.2 Probability of defibrillation success strongly depends on time when shock is applied

We now want to show how the success rate is time-dependent or, more specifically, that its time series shows rather high maxima which are not just small deviations around the mean success rate. We will then show that these peaks have a non-vanishing width. As already mentioned in the introduction, the aim of this study is to terminate arrhythmic cardiac excitation waves with perturbation parameters corresponding to low defibrillation field strengths, thus causing less pain and tissue damage due to the shocks, while keeping the termination success rate significantly higher than on average for the same energy. For the perturbations, we chose a fixed number of N_pert = 500 perturbation sites and a fixed amplitude of $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.5$ a.u. For this choice, the averaged DRC, Figure 3, yields a very low success rate of S(N_pert = 500) ≈ 4.3 % which corresponds to the success rate of applying shocks at random times (states) to the chaotic cardiac excitation. This is an unacceptably small value for clinical application.

Given the success rate distribution DRC, Figure 4, we analyze success rate time series S (t_appl) in terms of a whether they exhibit a significant peak structure. Compared to both DRCs shown in Figures 3, 4, we strongly increased the amount of realizations for the success rate analyses. First of all, we generated 100 realizations of the time evolution (i.e. 100 trajectories with different initial conditions), each being 10 s long, exhibiting chaotic spiral wave dynamics which were then analyzed time-wise in terms of their response to applied perturbations. The 100 chaotic trajectories were each sampled with a rate of 100.1 Hz, resulting in sequences of states. Each of these states was then perturbed individually 100 times using 100 different, randomly generated fixed spatial distributions of 500 perturbations sites (acting like virtual electrodes). In this way, a success rate was estimated for each of the 1,001 states (or points in time separated by dt = 10 ms) along the chaotic trajectory (CT) respresenting a fibrillation episode. The total amount of termination attempts realized is thus:

$\left(\text{no.\hspace{0.17em}of\hspace{0.17em}CT}\right)\times \left(\text{no.\hspace{0.17em}of\hspace{0.17em}}{t}_{\text{appl}}\right)\times \left(\text{no\hspace{0.17em}of\hspace{0.17em}conf.}\right)=100\times 1001\times 100\approx 10^7.$

First, we want to discuss the structure of the time series and the resulting implications. A typical example of such a time series of the success rate is given in Figure 5. As can be clearly seen, there is a well-established peak structure of the success rate and thus, a time-dependence of it. This means that there are short time intervals within which the termination of chaotic excitation waves is much more probable than on average. As a matter of fact, if external stimuli are applied to the most favorable state of the time series shown in Figure 5, which is at t_appl = 3.99 s, the probability to successfully terminate the unwanted dynamics is almost fifteen times higher than when averaging over all success rates of the 1,001 different states the trajectory consists of. On average, such a success rate would be reached only with twice as much perturbations, see Figure 3. The temporally averaged success rate for the trajectory shown in Figure 5 is $\bar{S}\approx 6.5$ % which agrees well with the order of value of the averaged DRC, Figure 3. The apparent difference (6.5 % vs. 4.3 %) is due to the fact that in the DRC, for each number of perturbation sites, 20 × 50 = 1,000 realizations are evaluated, whilst the time series shown in Figure 5 is composed of more than 100,000 termination attempts. Note that, for a sufficient number of realizations (perturbed states), the average termination success rate of a trajectory, $\bar{S}$ in Figure 5, and of randomly chosen single chaotic states, as for the averaged DRC, Figure 3, may be considered equivalent quantities and correspond to the termination success probability of defibrillating at random times.

FIGURE 5

FIGURE 5. Time series of the success rate, S(t_appl), for perturbing a trajectory of chaotic cardiac excitation. With a temporal resolution of dt_appl = 10 ms, the states along the trajectory are perturbed using 100 different spatial configurations of N_pert = 500 perturbation sites that were randomly chosen before starting the analysis and kept fixed along the trajectory. For each of the 100 defibrillation attempts, the success rate was evaluated according to Section 2 and the spiral wave dynamics termination success was estimated in terms of whether phase singularities still exist 500 ms after the perturbation. The temporal average of the success rate over the whole trajectory is $\bar{S}\approx 6.5\%$. The peaks have a certain width and are not just particularly favorable points in time for shocking, as shown exemplary for a higher resolved peak in Figure 6.

For all of the 100 trajectories, we found qualitatively similar success rate time series. Altogether, the main message of such peak structures for S(t_appl) is that, if perturbations are efficiently applied, i.e., at an appropriate time, one can yield success rates being, in parts significantly, higher than the temporal average, $\bar{S}$, over the trajectory. Another remarkable observation of the exemplary realization shown in Figure 5 (and others not shown here) is that the peaks possess a certain width. In order to illustrate this feature, we generated higher resolution plots of several of the peaks from the chaotic trajectories (i.e., again with N_pert = 500), using a time stepping of dt_appl = 0.1 ms (which requires the initial temporal integration step to be halved) for the termination attempts yielding a resolution being 10² times higher than used otherwise within the study. In Figure 6, we can exemplarily see a higher resolution of the largest peak of the time series indicated in Figure 5, i.e., at S(t_appl = 3.99 s) = 96 % which has a width of about 40 ms with respect to the mean value of the success rate $\bar{S}$. A similar width was also observed for other high-resolution success peaks (not shown here). The fact that high success rates do not occur within extremely narrow peaks but in relatively broad time windows is important with regard to their possible future use in practical experimental or clinical applications.

FIGURE 6

FIGURE 6. Higher-resolution plot of the peak highlighted in Figure 5 at 3.99 s for N_pert = 500, using smaller time steps for the application of local stimuli, dt_appl = 0.1 ms. The full width at half maximum (FWHM) is roughly 25 ms. However, the whole peak and thus, the time window with a success rate being (in parts significantly) higher than the average success rate $\bar{S}$ has a width of about 40 ms.

3.3 The peak structure is sufficiently robust to variations of perturbation parameters

Since the choices of the perturbation parameters, N_pert = 500 and $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.5$ a.u., are rather specific, in what follows, we will show that S(t_appl) is sufficiently robust to changes in these parameters, starting with an investigation of the robustness to changes in the number of local stimuli added to the system, N_pert. In order to do so, we exemplarily show a success rate time series of one out of the 100 trajectories perturbed with different N_pert ∈ {250, 450, 500, 550, 750, 1000}, respectively, evaluating again 100 different configurations of spatial distributions of perturbation sites for the same states, t_appl. In Figure 7, one can see the resulting time series of the success rate S(t_appl) if the very same trajectory is perturbed with these different numbers of perturbation sites. As N_pert is increased, the peak structure is, at least partially, preserved. Furthermore, most peaks that are preserved when applying more perturbations, grow in height. A good example where this can be observed is the highest peak for N_pert = 500 at around t_appl = 8.86 s. It does not only dominate the time series for N_pert = 500 but also the ones for the other numbers of perturbation sites used, growing in height the more N_pert are used. Only for the time series using 750 and 1,000 perturbation sites, the phenomenon of a few dominating peaks cannot be observed anymore because most of their peaks are generally very high. In terms of the hierarchy of peak heights, however, there does not seem to be a definite rule. The three peaks for N_pert = 250 appearing after 2 s have a different order in terms of their height if defibrillation attempts with more local stimuli are evaluated, for example. Finally, there are also examples where peaks present for a small number of perturbation sites seem to disappear as N_pert is increased. An example for this can be observed for the peak at 4 s for the time series with N_pert = 250 which disappears for N_pert ∈ {450, 500, 550} before finally recovering slowly as more perturbations sites are applied. In general, however, the corresponding temporal averages, $\bar{S}$, agree well with the averaged DRC in Figure 3, as it can be seen by taking a look at the monotonically increasing $\bar{S}$ as N_pert is increased.

FIGURE 7

FIGURE 7. Time series of the termination success rate of a single trajectory of chaotic cardiac excitation for different amounts of local stimuli applied N_pert to the system. $\bar{S}$ represents the corresponding average success rate of the time series S(t). Perturbations are applied as described in the caption of Figure 5.

In order to quantify these observations, we calculated the Spearman and Pearson correlation coefficients, r_s and ρ_p, respectively, between all combinations of time series given in Figure 7. The Spearman correlation is a measure for how well a monotonic function can approximate the relationship between two samples or, in this case, between two time series. In particular, the Spearman coefficients between S(t_appl) for N_pert = 500 and S(t_appl) of the closest numbers investigated here, i.e., 450 and 550, are 0.78 and 0.85, respectively. The complete set of Spearman and Pearson correlation coefficients between the success time series for different numbers of perturbation sites is given in the Supplementary Material. The Pearson correlation coefficient measures the linear dependence of two variables. For small changes in N_pert = 500, so again N_pert = 450 and N_pert = 550, the Pearson correlation coefficient of the corresponding time series S(t_appl) is even higher, ρ_p = 0.93 in both cases, respectively. The relatively large value of the Spearman coefficient and the even larger value for the Pearson coefficient between time series perturbed with similar numbers of perturbation sites, agrees well with the qualitative observation that, most of the peaks are preserved with growing heights, whereas some of them seem to vanish as N_pert is increased. Hence, we may conclude that the peak structure is sufficiently robust, especially to small variations in N_pert.

We now want to show a similar robustness of S(t_appl) to changes in the second perturbation parameter, the amplitude of local stimuli, $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}$. Therefore, we repeated the previous investigation, this time by varying the magnitude of the perturbations and not their number, using $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}\in \{0.15,0.3,0.45,0.5,0.55,0.7\}$ a.u. In principle, the situation is similar to the time series with varying N_pert, i.e., rather large parts of the peak structure is preserved if the perturbation amplitude is increased, as illustrated in Figure 8. However, for very small amplitudes, $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.15$ a.u., this does not apply. Taking a look at the peaks of this success rate time series and compare it with the peaks of S(t_appl) with a perturbation magnitude twice as large, we see that most of the peaks are attenuated which is counter-intuitive to the observations made so far. When taking a look at the correlation coefficients between the S(t_appl) perturbed with the magnitude used here, $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.5$ a.u., with its two neighboring magnitude values investigated, $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.45$ a.u. and $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.55$ a.u., however, we obtain values similar to the variation of N_pert around 500 perturbations. The Spearman correlation coefficients are r_s = 0.8 and r_s = 0.87, respectively, whereas the Pearson correlation coefficient is the same for varying $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}=0.5$ a.u. in both directions, namely ρ_p = 0.95. The complete set of correlation coefficients between all trajectories with varying perturbation magnitude is given in the Supplementary Material as well. This demonstrates that also variations of the perturbation amplitude chosen for this study are sufficiently robust in terms of the success rate time series. The fact that the peak structure is not completely preserved when changing N_pert and $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}$ is very likely a consequence of the chaotic nature of the trajectory and hence, probably inevitable. However, the general finding here is that the peak structure is relatively robust to changes in the perturbation parameters.

FIGURE 8

FIGURE 8. Time series of the termination success rate of a single chaotic trajectory for different magnitudes of the perturbations, $\mathrm{\Delta }{V}_{\text{m}}^{2\times 2}$, applied to the system. $\bar{S}$ represents the average success rate of the respective time series.

3.4 Relation between peak height and waiting times

Even though most trajectories do exhibit states with relatively large success rates, these peaks do not occur very frequently, as illustrated by the time series shown so far. We therefore show here how waiting times for peaks increase dramatically the larger we require a peak to be and thus, the success rate for defibrillation. For each state of a trajectory investigated, given by its time of occurrence tⁱ, the waiting time is defined as the distance in time from this state to its closest temporally forward peak of the required height h

$t_{\text{wait}}^i=t_h^i-t^i.(10)$

Note that, if tⁱ itself meets the condition (i.e., is the point in time at which a peak of given height h occurs), then $t_{\text{wait}}^i=0$ ms. Furthermore, states after the last $t_h^i$ of a specific trajectory are not considered as there is no subsequent closest temporally forward peak (state) for them to refer to. In Figure 9, the distribution of the waiting times of each sampled state the 100 trajectories are composed of is plotted semi-logarithmically for h ∈ {20, 40, 60} %. The waiting times for a peak with a termination success rate of 20 % are usually not longer than 2 s and the fraction of waiting times smaller than 1 s is 0.94. For 60 %, however, this fraction reduces dramatically down to 0.26. In the context of defibrillation, this means that one usually does not have to wait long for a time interval inside which the success rate is multiple times larger than on average, but the success rates are the lower the less one wants to wait for the defibrillation to be applied, on average.

FIGURE 9

FIGURE 9. Probability density of waiting times for the next peak (with specific minimal heights, respectively) to occur inside one trajectory, averaged over all states of 100 trajectories that do not correspond to a peak of the success rate, see Eq. 10.

Finally, we show here how the peak structure, S(t_appl), is characterized in terms of its frequency by comparing the frequency spectra obtained via Fourier transformation (FFT) of both, the success rate and the mean of the membrane potential time series, S(t_appl) and ${\bar{V}}_{\text{m}}(t)$ where, for simplicity, t_appl = t. In order to do so, we normalized each of the 100 time series, S(t) and ${\bar{V}}_{\text{m}}(t)$, with respect to their corresponding maximal value, denoted by $\widehat{S}(t)$ and ${\stackrel{^̄}{V}}_{\text{m}}(t)$ (not to be confused with the hat notation for Fourier transforms). We then calculated the FFT of these time series, evaluated their modulus and finally averaged over the moduli of the FFTs associated with each trajectory, yielding the amplitude spectra of both quantities. The comparison of these frequency spectra of ${\bar{V}}_{\text{m}}(t)$ and S(t) is shown in Figure 10.

FIGURE 10

FIGURE 10. Mean amplitude spectra of the success rate and the membrane voltage averaged over 100 realizations (trajectories). Note that the hat denotes the normalization of both quantities before the carrying out the FFT in order to obtain comparable magnitudes and does not itself denote an FFT.

The strongest frequency components of the 100 time series ${\bar{V}}_{\text{m}}(t)$ occur between 2 and 7 Hz. This shape of the spectrum makes it difficult to specify a single frequency as “the” dominant frequency, but the range [2 Hz, 7 Hz] agrees with experimental data on ventricular fibrillation in human hearts, for example in Pandit and Jalife (2013). If a dominant frequency can be clearly identified, its inverse corresponds to the (average) rotation period of the spiral waves Mandapati et al. (1998) and thus provides a relevant time scale characterizing the chaotic wave dynamics. An interesting observation is that the peak of the success rate spectrum at 7.5 Hz occurs at the upper limit of the spectrum of the membrane voltage. This means that within each rotational period of even the fastest spiral waves (on average), there is a state for which the application of perturbation would be particularly successful in terminating the spatiotemporal chaos.

4 Discussion and conclusion

Defibrillation shocks used to terminate the complex electrical excitation dynamics in the heart during ventricular fibrillation must be of high energy to achieve a high success rate in terminating the arrhythmia, resulting in pain and damage to heart tissue Holden et al. (2019). In practice, these shocks are applied at random times during the fibrillation. Therefore, in this study, we addressed the question of whether there are more favorable times (or states) in the temporal evolution of an arrhythmia at which even medium- or low-energy shocks can terminate the chaotic state corresponding to fibrillation with relatively high probablility. To investigate this question, we performed extensive numerical simulations using the Fenton-Karma model and an implementation of multi-site pacing by virtual electrodes. With these simulations, we were not only able to reproduce the characteristic sigmoidal shape of the dose-response curve, but also found strong fluctuations of the termination success rate during the time evolution of the complex dynamics of cardiac activity resulting in well-established peak structures of success rate time series. In these simulations, it was found that even relatively low energy shocks could terminate the chaotic dynamics if they were applied in short time windows in which the termination success was particularly high. If such fluctuations of termination success could be predicted, they might be exploited in future defibrillation protocols for terminating arrhythmias with low or medium pulse energies.

Since the underlying model represents a simple model for cardiac excitation dynamics, the simulations made here to study success rate time series of trajectories representing arrhythmias only qualitatively describe the dynamics of a real heart. First of all, the Fenton-Karma model is of low to moderate complexity and thus, may not incorporate cardiac dynamics with all desired details. Furthermore, out of several parameter sets suggested for the Fenton-Karma model, we only used a specific one. Thus, the results found in the study presented here need to be confirmed in future investigations. As a starting point, one could repeat the simulations using different parameter sets of the Fenton-Karma model Fenton et al. (2002), corresponding to different properties of cardiac dynamics. Besides being computationally more expensive, the results presented should also be verified using more detailed (ionic) models of cardiac excitation, such as the 4-variable Bueno-Orovio-Cherry-Fenton model Bueno-Orovio et al. (2008), the 8-variable Beeler-Reuter model Beeler and Reuter (1977) or the 17-variable Ten Tusscher-Noble-Noble-Panfilov model ten Tusscher et al. (2004). Future studies may also consider more realistic geometries of the (human) heart Fenton et al. (2005) including fiber orientation Doste et al. (2019) influencing the excitation propagation. Heterogeneities, such as blood vessels and collagen, may not only act as virtual electrodes but also have an impact on the spatio-temporal dynamics Bittihn et al. (2017); Rappel et al. (2022). Modelling this impact in detail requires the use of bidomain models differentiating between intra- and extracellular space. For a more realistic modeling of virtual electrodes, one may use stimuli based on local current injection Feng et al. (2022) and consider perturbation sites of different sizes and shapes where different excitation thresholds due to the strength-extent relation have to be taken into account Bezekci et al. (2015). And again, in all these simulations it has to be checked whether the peaks of high success rate are (still) broad enough to practically apply external shocks precisely within these time frames. Finally, if the findings of the presented study can be confirmed in more detailed simulations of arrhythmic cardiac dynamics and fibrillation, the next challenge will be to predict the peaks in the success rate from measurable observables of the heart.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

JS, UP, TL, and SL conceived the study, JS performed the simulations and analyzed the data. JS, UP, TL, and SL wrote the manuscript.

Funding

The authors acknowledge support by the Max Planck Society (MPG), the German Centre for Cardiovascular Research (DZHK) e.V. and the European High-Performance Computing Joint Undertaking EuroHPC under grant agreement No 955495 (MICROCARD).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2022.1007585/full#supplementary-material

References

Beeler, G. W., and Reuter, H. (1977). Reconstruction of the action potential of ventricular myocardial fibres. J. Physiology 268, 177–210. doi:10.1113/jphysiol.1977.sp011853

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezekci, B., Idris, I., Simitev, R. D., and Biktashev, V. N. (2015). Semianalytical approach to criteria for ignition of excitation waves. Phys. Rev. E 92, 042917. doi:10.1103/PhysRevE.92.042917

PubMed Abstract | CrossRef Full Text | Google Scholar

Bittihn, P., Berg, S., Parlitz, U., and Luther, S. (2017). Emergent dynamics of spatio-temporal chaos in a heterogeneous excitable medium. Chaos 27, 093931. doi:10.1063/1.4999604

PubMed Abstract | CrossRef Full Text | Google Scholar

Bittihn, P. (2015). Complex structure and dynamics of the heart. 1 edn. Switzerland: Springer International Publishing.

Google Scholar

Bittihn, P., Luther, G., Bodenschatz, E., Krinsky, V., Parlitz, U., and Luther, S. (2008). Far field pacing supersedes anti-tachycardia pacing in a generic model of excitable media. New J. Phys. 10, 103012. doi:10.1088/1367-2630/10/10/103012

CrossRef Full Text | Google Scholar

Boegli, M., Blanc, O., Virag, N., Vesin, J.-M., and Kappenberger, L. (2000). “Study of the defibrillation process in a computer model of human atria,” in Proceedings of the 22nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Cat. No.00CH37143), Chicago, IL, USA (IEEE), 1848–1849. doi:10.1109/IEMBS.2000.900446

CrossRef Full Text | Google Scholar

Bragard, J., Elorza, J., Cherry, E., and Fenton, F. (2013). “Validation of a computational model of cardiac defibrillation,” in Computing in Cardiology 2013, Zaragoza, Spain (IEEE), 851–854.

Google Scholar

Bueno-Orovio, A., Cherry, E., and Fenton, F. (2008). Minimal model for human ventricular action potentials in tissue. J. Theor. Biol. 253, 544–560. doi:10.1016/j.jtbi.2008.03.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Cates, A. W., Wolf, P. D., Hillsley, R. E., Souza, J. J., Smith, W. M., and Ideker, R. E. (1994). The probability of defibrillation success and the incidence of postshock arrhythmia as a function of shock strength. Pacing Clin. Electrophysiol. 17, 1208–1217. doi:10.1111/j.1540-8159.1994.tb01487.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chugh, S. S. (2017). Sudden cardiac death in 2017: Spotlight on prediction and prevention. Int. J. Cardiol. 237, 2–5. doi:10.1016/j.ijcard.2017.03.086

PubMed Abstract | CrossRef Full Text | Google Scholar

Datseris, G., and Parlitz, U. (2022). Nonlinear dynamics - a concise introduction interlaced with code. Springer.

Google Scholar

Davy, J.-M., Fain, E. S., Dorian, P., and Winkle, R. A. (1987). The relationship between successful defibrillation and delivered energy in open-chest dogs: Reappraisal of the “defibrillation threshold” concept. Am. Heart J. 113, 77–84. doi:10.1016/0002-8703(87)90012-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Doste, R., Soto-Iglesias, D., Bernardino, G., Alcaine, A., Sebastian, R., Giffard-Roisin, S., et al. (2019). A rule-based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. Int. J. Numer. Method. Biomed. Eng. 35, e3185. doi:10.1002/cnm.3185

PubMed Abstract | CrossRef Full Text | Google Scholar

Fast, V. G., and Cheek, E. R. (2002). Optical mapping of arrhythmias induced by strong electrical shocks in myocyte cultures. Circ. Res. 90, 664–670. doi:10.1161/01.RES.0000013403.24495.CC

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, X., Yin, X., Wen, J., Wu, H., and Gao, X. (2022). Removal of spiral turbulence by virtual electrodes through the use of a circularly polarized electric field. Chaos 32, 093145. doi:10.1063/5.0102031

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F., Cherry, E., Hastings, H., and Evans, S. (2002). Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos (Woodbury, N.Y.) 12, 852–892. doi:10.1063/1.1504242

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F. H., Cherry, E. M., Karma, A., and Rappel, W.-J. (2005). Modeling wave propagation in realistic heart geometries using the phase-field method. Chaos 15, 013502. doi:10.1063/1.1840311

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F., and Karma, A. (1998). Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos (Woodbury, N.Y.) 8, 20–47. doi:10.1063/1.166311

PubMed Abstract | CrossRef Full Text | Google Scholar

Holden, A. V., Begg, G. A., Bounford, K., Stegemann, B., and Tayebjee, M. H. (2019). Phase entrainment of induced ventricular fibrillation: A human feasibility and proof of concept study. J. Atr. Fibrillation 12, 2217. doi:10.4022/jafib.2217

PubMed Abstract | CrossRef Full Text | Google Scholar

Iyer, A., and Gray, R. (2001). An experimentalist’s approach to accurate localization of phase singularities during reentry. Ann. Biomed. Eng. 29, 47–59. doi:10.1114/1.1335538

PubMed Abstract | CrossRef Full Text | Google Scholar

Janardhan, A. H., Gutbrod, S. R., Li, W., Lang, D., Schuessler, R. B., and Efimov, I. R. (2014). Multistage electrotherapy delivered through chronically-implanted leads terminates atrial fibrillation with lower energy than a single biphasic shock. J. Am. Coll. Cardiol. 63, 40–48. doi:10.1016/j.jacc.2013.07.098

PubMed Abstract | CrossRef Full Text | Google Scholar

Koplan, B., and Stevenson, W. (2009). Ventricular tachycardia and sudden cardiac death. Mayo Clin. Proc. 84, 289–297. doi:10.1016/S0025-6196(11)61149-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwaku, K. F., and Dillon, S. M. (1996). Shock-Induced depolarization of refractory myocardium prevents wave-front propagation in defibrillation. Circ. Res. 79, 957–973. doi:10.1161/01.RES.79.5.957

PubMed Abstract | CrossRef Full Text | Google Scholar

Lilienkamp, T. (2018). A state space odyssey — the multiplex dynamics of cardiac arrhythmias. Göttingen: University of Göttingen. Dissertation.

Google Scholar

Lilienkamp, T., and Parlitz, U. (2020). Terminating transient chaos in spatially extended systems. Chaos 30, 051108. doi:10.1063/5.0011506

PubMed Abstract | CrossRef Full Text | Google Scholar

Luther, S., Fenton, F., Kornreich, B., Squires, A., Bittihn, P., Hornung, D., et al. (2011). Low-energy control of electrical turbulence in the heart. Nature 475, 235–239. doi:10.1038/nature10216

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandapati, R., Asano, Y., Baxter, W. T., Gray, R., Davidenko, J., and Jalife, J. (1998). Quantification of effects of global ischemia on dynamics of ventricular fibrillation in isolated rabbit heart. Circulation 98, 1688–1696. doi:10.1161/01.CIR.98.16.1688

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandit, S. V., and Jalife, J. (2013). Rotors and the dynamics of cardiac fibrillation. Circ. Res. 112, 849–862. doi:10.1161/CIRCRESAHA.111.300158

PubMed Abstract | CrossRef Full Text | Google Scholar

Press, W. H., Teukolsky, A. S., Vetterling, T. W., and Flannery, B. P. (2007). Numerical recipes 3rd edition: The art of scientific computing. USA: Cambridge University Press.

Google Scholar

Pumir, A., and Krinsky, V. (1999). Unpinning of a rotating wave in cardiac muscle by an electric field. J. Theor. Biol. 199, 311–319. doi:10.1006/jtbi.1999.0957

PubMed Abstract | CrossRef Full Text | Google Scholar

Rappel, W.-J., Krummen, D. E., Baykaner, T., Zaman, J., Donsky, A., Swarup, V., et al. (2022). Stochastic termination of spiral wave dynamics in cardiac tissue. Front. Netw. Physiol. 2, 809532. doi:10.3389/fnetp.2022.809532

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H. W. J., Noble, D., Noble, P. J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. Heart Circ. Physiol. 286, H1573–H1589. doi:10.1152/ajpheart.00794.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N. A., Skouibine, K., and Aguel, F. (1998). The role of cardiac tissue structure in defibrillation. Chaos 8, 221–233. doi:10.1063/1.166299

PubMed Abstract | CrossRef Full Text | Google Scholar