This section provides the mathematical framework to describe the myocardium, which consists of a network of electrically and mechanically coupled cardiac cells.

Apart from defibrillation attempts with a single pulse (SinglePulse), numerous low energy pacing protocols exist that are based on the application of pulse sequences derived from variables of the frequency spectrum of the underlying dynamics. For example, pulses can be applied at a constant frequency (EquiDist) (Fenton et al., 2009; Luther et al., 2011; Janardhan et al., 2014; Lilienkamp et al., 2022a) or at decelerating pulse frequencies (ADP) (Lilienkamp et al., 2022b). Both protocols share the characteristic that all pulse frequencies (and thus pulse timings) are calculated based on observations of the system in the last seconds prior to the first pulse. During the pacing process, the dynamics of the system is neglected, as the predetermined pulse frequencies are kept fix. However, Steyer et al. showed that the exact timing of a single pulse plays a major role for a successful termination of a chaotic state (Steyer et al., 2023). Additionally, the accuracy of the calculated pulse frequencies is inherently limited by the time period for which the system is observed prior to pacing (B Traykov et al., 2012), as ADP and EquiDist are based on the frequency spectrum of the underlying dynamics and, respectively, its dominant frequency. Furthermore, Buran et al., 2023 observed that successful low energy defibrillation protocols are characterized by a coordinated interplay between the applied pulses and other observables of the underlying dynamics, such as the refractory boundary length, and proposed a feedback-controlled strategy based on this observable (Buran, 2023).

We intend to contribute to the advancement of feedback-controlled protocols by utilizing the mean value of the transmembrane potential V _m V m ̄ = 1 N 2 ∑ i = 1 N ∑ j = 1 N V m i,j (4) as the observable to calculate pulse timings. Figure 2A illustrates, that local minima of V m ̄ are followed by local maxima of v ̄ and w ̄ for the FK model. Intuitively, these local maxima of the mean values of the gating variables indicate times at which the system exhibits a high degree of excitability. Pulses applied to the system at these points in time are therefore more likely to excite a larger fraction of the domain than pulses at other times. The same phenomenon can be observed in the AP, BOCF, and MSA model (see Supplementary Material). This observation motivates our choice to use local minima of V m ̄ as marker points to define the timings, at which pulses should be applied. We define the upward deflection time Δt _up as the time between a local minimum and the consecutive local maximum in V m ̄ (compare Figure 2B). The median Δ t median = Δ t up ̃ of all upward deflection times over all initial conditions of a specific model of cardiac dynamics is calculated based on the last 2 seconds of V m ̄ prior to the first pulse. Therefore, we determined the median upward deflection time for each model of cardiac dynamics, separately.

In Local Minima Pacing (LMP), a pulse is applied shortly after a local minimum of V m ̄ is detected. Subsequently, the system evolves autonomously, the feedback of the system to the given pulse is measured and a consecutive pulse is applied after the next minimum is detected. The key parameter of LMP is the temporal distance between a local minimum and the application of a pulse, which we define with respect to the median upward deflection time. A pulse is applied if ∂ t V m ̄ ( t ) > 0 , ∀t ∈ [t _c − a Δt _median, t _c] and t _c − t _last,p > t _pause, where t _c is the current time, a acts as a scaling factor, t _last,p is the time of the last LMP pulse, and t _pause prevents consecutive pulses during an upward deflection of V m ̄ initiated by a LMP pulse. Viable parameter selections for a are analysed in Section 3.2. LMP stops if a maximum number of pulses max _pulses is reached. By definition, however, it also stops if max i,j V m i,j < ϵ and therefore if all spiral wave patterns are terminated successfully. A description of the LMP approach is also given in Algorithm 1.

Algorithm 1

Local Minima Pacing - Pseudocode.

Input: V m ̄ Mean value of the transmembrane potential

In this section, we begin to systematically study LMP, using a specific parameterization of it to understand how LMP pulses influence the dynamics of a system that exhibits spiral wave patterns. Figure 3 demonstrates the timing and impact of LMP pulses on the membrane potential (see Figure 3A) and the mean value of the membrane potential (Figure 3B) of the same state for the FK model in comparison to a temporal evolution without perturbations. For the analysis of the LMP algorithm presented in this subsection, we set the scaling factor for the pulse timings to a = 0.2, implying that pulses were applied at 20% of the median upward deflection time Δt _median after a local minima in V m ̄ was detected.

To compare the efficiency of LMP with other previously published pacing protocols, namely, conventional defibrillation (SinglePulse), Equidistant Pacing (EquiDist), and Adaptive Deceleration Pacing (ADP), we calculated dose-response curves, shown in Figure 4, for the AP model (see Figure 4A), the BOCF model (see Figure 4B), the FK model (see Figure 4C), and the MSA model (see Figure 4D), respectively. This way, a relation between the pulse amplitude A _max and the success rate to terminate spiral wave dynamics can be determined. The dominant frequency and the frequency spectrum used in EquiDist and ADP were determined on 5 s time series of the underlying dynamics, as described by Lilienkamp et al., 2022b. The pacing frequency of EquiDist(c) was set to c f _dom, with c > 0 as a scaling factor. The local minimum of the EquiDist (0.8) dose response curve of the FK model at A _max ≈ 1.2 [a.u.] (see Figure 4C), was previously discussed by Lilienkamp et al., 2022a and we can confirm the reasoning of Buran et al., 2017 that later pulses may reinitiate fibrillation after successful termination.

We investigated how the scaling parameter a and therefore the temporal distance between a local minimum and a LMP pulse affects the overall performance of the LMP protocol. The tests we carried out cover a range of 10%–100% of the median upward deflection time Δt _median as the timing of the LMP pulses. Motivated by decelerating pulse frequencies in the ADP protocol, we allow for increasing or decreasing scaling factors in LMP and incorporate an adjustment of the scaling factor a. After a LMP pulse is applied, the scaling factor is adjusted using a = clamp (a b, Δt _start/Δt _media, Δt _end/Δt _median), with b > 0 as an adjustment factor for the pulse timings during LMP. Here, Δt _start describes the temporal difference between the first and second LMP pulse and Δt _end between the last two pulses, respectively. This update rule is independent of the maximum number of pulses max _pulses, permits an adaptive temporal distance between a local minimum in V m ̄ and a LMP pulse but retains the central concept of LMP: an excitation of an increasing fraction of the domain is more likely right after local minima in V m ̄ .

Figures 5A–D show a quantitative comparison of the pulse timings with respect to the required pulse amplitudes A _max to obtain a 90% success rate. To account for the non-linearity of the update rule, we set b = 0.8 if Δt _start > Δt _end and b = 1.25 if Δt _start < Δt _end, where Δt _start denotes the temporal difference between the first local minimum and the first LMP pulse in ms, Δt _end for the last minimum and LMP pulse, respectively. We implemented ADP and EquiDist to perform pulse sequences of ten pulses. Therefore, the maximum number of LMP pulses was set to max _pulses = 10. Figures 5A–D show, that for each model multiple combinations of Δt _start and Δt _end exist, for which the required E90-energy of LMP is lower compared to all other protocols (marked with purple crosses). However, the key point of these plots is, that LMP with the combination of Δt _start = Δt _end = 0.2 Δt _median offers an energy reduction for all models of cardiac dynamics described in Section 2.1.2. To achieve a 90% success rate, we observed a reduction in the required amplitude A _max of 28% for the AP model, 0.7% for the BOCF model, 56% for the FK model, and 18% for the MSA model compared to the ADP protocol in our numerical simulations.

Previously published pacing protocols like EquiDist and ADP are mainly parameterized through the pulse timings of a pulse sequence (the temporal distance between consecutive pulses as depicted in Figure 3B). The proposed LMP approach, however, is parameterized by the temporal difference between a local minimum in V m ̄ and a LMP pulse. To deepen our insight into the LMP approach we fixed the scaling factor to a = 0.2, performed an a posteriori analysis of the temporal distances (periods) between LMP pulses at the respective E₉₀-amplitude and compared them with pulse timings from the EquiDist and ADP protocols. The relative occurrences for a given temporal distance are shown in Figure 6, where models of cardiac dynamics are displayed in separate rows and different pacing protocols in columns, respectively. The dominant period and the underlying frequency spectra of the ten initial conditions we simulated for each model of cardiac dynamics vary slightly. Therefore, the pacing periods derived by EquiDist and ADP show bundled patterns (see first and second column of Figure 6). LMP was parameterized on averages of Δt _median of each model, and not specified on observations of each initial condition prior to pacing like EquiDist and ADP. On average, the resulting LMP periods show a decelerating pattern (compare second and third column of Figure 6), because the action potential durations increase with a growing number of LMP pulses (compare Figure 3B) but saturate at a level, that represents the action potential duration of a plane wave. It should be noted that LMP is not limited to strictly decelerating sequences of pulses. We also observe single pulse sequences that show monotonously increasing pulse periods and sharply decreasing periods in the last two pulses. Furthermore, LMP periods cover a broader range in each pulse pair compared to ADP. This observation highlights the adaptive nature of LMP, which incorporates the feedback of the underlying dynamics during pacing.

Distinct gaps in the period distribution in single pulse pairs, as shown in Figure 6I between 170–180 ms, indicate occurrences of saddle points instead of local minima of V m ̄ in a subset of termination attempts. Filled period bins at 0 ms for given pulse pairs indicate successful termination without applying further pulses. For the AP (see Figure 6C), the BOCF (see Figure 6F), the FK (see Figure 6I), and the MSA (see Figure 6L) model, 50% of all successful termination attempts only required nine, five, eight, and seven pulses, respectively. Therefore, LMP not only permits the usage of lower amplitudes but also potentially reduces the number of pulses required for successful termination of the underlying dynamics. This is not the case for previously discussed termination strategies like sequences of equidistant pulses or the ADP protocol. Buran et al., 2017 describe a reinitialisation of spiral waves by later pulses in EquiDist after successful termination, which we can confirm for EquiDist and ADP. Due to its design, LMP is not susceptible to behaviour.

The feedback-controlled pacing protocol proposed by (Buran, 2023) determines pulse timings based on minima in the length of the refractory boundary L _RB. These minima, however, do not always coincide with the minima in V m ̄ . This observation is in line with the findings of (Buran, 2023): The correlation between the length of the refractory boundary L _RB and the excitable fraction of the domain either do not exist at all or only in a weak manner, for systems exhibiting spatiotemporal chaos and unstable spirals.

As described in Section 2.1.4, the coverage area A _cov describes the relative fraction of the simulation domain that is covered by AVEs and was modelled very differently in previous publications (Gray, 2002; Buran et al., 2017; Lilienkamp et al., 2022b). Like other modelling parameters, the coverage area A _cov can only be approximated. There is no guarantee that this approximation applies to every individual, which is why it is important to perform a sensitivity analysis in order to be able to assess the statements based on the underlying experiments (Pathmanathan et al., 2019; Salvador et al., 2023). Therefore, we extended our investigations, presented in the previous subsections, and performed termination attempts at 80% and 60% of the coverage rates stated in Table 1 and examined how sensitive dose response curves of different pacing protocols and model of cardiac dynamics depend on different values of A _cov. Figure 7 illustrates the effect of lower coverage areas A _cov and therefore lower numbers of AVEs N _ave, as we kept the size of single AVEs constant as stated in Table 1. Lower coverage areas require higher amplitudes A _max to achieve similar success rates. Figures 7A–D emphasise that the aforementioned results change quantitatively, but not qualitatively, with lower coverage rates A _cov. Figures that include dose response curves for EquiDist and ADP are given in the Supplementary Material.

The numerical experiments presented in this work were performed on a two dimensional grid, on simplified models of cardiac dynamics, and without the incorporation of the muscle fibre orientation. Although the proposed approach is not yet designed to target a specific type of arrhythmia, by choosing a two-dimensional grid, our conclusions can be more directly applied to the thin, quasi-two-dimensional atrial wall (Luther et al., 2011). In particular, with respect to the approximation of the mean value of the transmembrane potential, further studies are needed to evaluate the performance of the proposed approach and the approximation of this observable on the surface of a three-dimensional grid representing the thick ventricular wall.

Furthermore, we neglected the coupling between cardiac mechanics and electrical dynamics (Ambrosi et al., 2011) in the myocardium, which influences the stability of spiral/scroll waves (Hu et al., 2013), and the shape of the action potential, especially during pacing (Hazim et al., 2021). We refer to successful defibrillation as the absence of chaotic patterns in the transmembrane potential, however, mechanical synchrony and function of the heart do not necessarily require electrical synchrony (Leclercq et al., 2002). In a further study, the proposed LMP approach should be tested in this regard.

The effect of unpinning stable spiral waves from tissue heterogeneities as a mechanism to terminate cardiac arrhythmias (Ripplinger et al. (2006)), has been studied extensively in numerical (Pumir and Krinsky 1999) and experimental (Gray and Chattipakorn, 2005; Ripplinger et al., 2006; Shajahan et al., 2016) setups.

We assume homogeneous and isotropic diffusion without conduction heterogeneities of the given scale. We therefore assume absence of scar tissue and neglect the influence of the fibre orientation (Hooks et al. 2007) in the myocardium. The interplay between LMP with the effect of unpinning, scar tissue and the influence of fibre orientation needs to be further investigated. In a future study addressing these limitations, it will be necessary to base the simulations on experimental data from specific mammals or humans.

Heterogeneities in the electrical conductivity of the myocardium are responsible for the depolarization of the transmembrane potential induced by an external electric field (Efimov et al., 2000; Pumir et al., 2007; Luther et al., 2011; Bittihn et al., 2012). In this study, this underlying mechanism is modelled by the application of local current injections at predefined locations within the simulation domain (Lilienkamp et al., 2022b), referred to as artificial virtual electrodes (AVEs). By employing the concept of AVEs, we accept that in particular the first milliseconds of the interaction between an external electric field and the myocardium is simplified. In previous studies, the utilization of AVEs proved effective in replicating the amplitude reductions of applying sequences of pulses at constant pulse timings in contrast to conventional defibrillation (Lilienkamp et al., 2022a; Lilienkamp et al., 2022b), which was also shown in other numerical (Buran et al., 2017) and experimental (Fenton et al., 2009) studies. Furthermore, applying sequences of low energy pulses may lead to the presence of local minima in the dose response curves for different models of cardiac dynamics (Buran et al., 2017), which can also be reproduced using the concept of AVEs (Lilienkamp et al., 2022a; Lilienkamp et al., 2022b). The aim of this study is therefore not to predict the exact shock strengths required by different pacing protocols, but rather to make qualitative statements about possible amplitude reductions of the protocols with respect to each other. However, this concept enables a grid spacing (Δx, N, and Δt) to perform a large number of simulations and thus a statistical analysis of the presented methods on multiple models of cardiac dynamics.

In this study, we investigated a novel approach as a promising alternative for low-energy defibrillation of cardiac arrhythmias like ventricular fibrillation. Multiple approaches, like equidistant pacing (Fenton et al., 2009; Luther et al., 2011) and adaptive deceleration pacing (Lilienkamp et al., 2022b), have been proposed to reduce the energy required for successful defibrillation. However, an energy reduction that enables defibrillation without side effects remains an open challenge.

In this study, we showed that pulses shortly after local minima in the mean value of the transmembrane potential can provide a significant energy reduction compared to recently published pacing protocols. Although the dynamics of the underlying system during pacing is incorporated in this approach, a single parameterization (scaling factor a = 0.2) of the proposed algorithm could be identified that provides a reduction of the energy for all investigated models of cardiac dynamics. This observation suggests that this parameterization may also result in energy reductions in future ex vivo experiments. Moreover, we provide a sensitivity analysis to demonstrate stability of the drawn conclusions with respect to a modelling parameter that was not selected uniformly in earlier publications. Compared to EquiDist and ADP, LMP potentially reduces the number of pulses and therefore the total electrical energy applied required to successfully terminate the underlying dynamics. Furthermore, as LMP incorporates the feedback of the underlying dynamics, a broader frequency band within each pulse pair is covered. This observation highlights the adaptive nature of this approach with respect to the underlying dynamics during pacing.

LMP solely relies on the measurement of the mean value of the transmembrane potential, an observable that can be measured very easily in numerical simulations, compared to the length of the refractory boundary (used in the feedback-controlled pacing algorithm proposed by Buran, 2023). As the measurement or approximation of both observables remains an open challenge in ex vivo experiments, we suggest testing the feasibility and reliability of these approximations in future experiments. In order to transfer this approach to more realistic numerical or even living-heart experiments, the mean value of the transmembrane potential needs to be computed, which remains an open challenge. We propose to reconstruct and approximate this observable by reconstructing the electrical activity of the heart from densely sampled body-surface potentials (Bear et al., 2018; Xie and Yao, 2022) in future studies. Of particular interest will be the noise induced by these approximation techniques and its influence on the minima and maxima of the mean transmembrane potential. The application of controlled electrical stimuli during LMP may further complicate the live measurement of the mean value of the transmembrane potential, e.g., by interactions between the external electric field and the ECG electrodes.

Characteristic features of atrial (AF) and ventricular fibrillation (VF) cover a broad range, depending on the biological species, medical condition and other factors. The mean dominant frequency of the underlying dynamics range from f _dom = 6.2 ± 0.2 Hz in the absence of cardiogenic shock and f _dom = 4.0 ± 0.2 Hz in the presence of cardiogenic shock for human VF (Stewart et al., 1992), f _dom = 5.7 ± 1.2 Hz for human AF (Bollmann et al., 1998), f _dom = 9.4 ± 2.6 Hz for AF in isolated sheep hearts (Skanes et al., 1998) to f _dom = 13.2 ± 3.4 Hz in isolated rabbit hearts (Chen et al., 2000). With regard to the mean dominant frequencies, this should provide a motivation for our choice of parameterizations of the models investigated, as these cover a range between f _dom = 2.65 Hz (BOCF) and f _dom = 22.46 Hz (AP). Furthermore, the mean number of phase singularities in our simulations range from N _PS = 10.74 (FK) to N _PS = 22.11 (AP). Although N _PS = 22.11 phase singularities clearly only applies to outliers, this range corresponds to recent numerical and experimental studies on pig and human hearts, respectively (Nash et al., 2006; Ten Tusscher et al., 2007).

Therefore our findings are not tailored to a specific model of cardiac dynamics, a parameterization of it, a single biological species, or one specific type of cardiac arrythmia but rather show a basic property inherent in all of these models, that can be utilized for defibrillation. Therefore, the exact energy reduction achievable with this approach, particularly taking into account the species, medical condition and type of cardiac arrhythmia, is therefore to be determined in future studies. The results presented, however, strongly indicate that LMP is a viable option to significantly lower the negative side effects associated with conventional defibrillation techniques.