^{1}Computational Physics for Life Science, Nuremberg Institute of Technology Georg Simon Ohm, Nuremberg, Germany^{2}Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany^{3}German Center for Cardiovascular Research (DZHK), Partner Site Göttingen, Göttingen, Germany^{4}Institute for the Dynamics of Complex Systems, Georg-August-Universität Göttingen, Göttingen, Germany^{5}Institute of Pharmacology and Toxicology, University Medical Center Göttingen, Göttingen, Germany

Current treatments of cardiac arrhythmias like ventricular fibrillation involve the application of a high-energy electric shock, that induces significant electrical currents in the myocardium and therefore involves severe side effects like possible tissue damage and post-traumatic stress. Using numerical simulations on four different models of 2D excitable media, this study demonstrates that low energy pulses applied shortly after local minima in the mean value of the transmembrane potential provide high success rates. We evaluate the performance of this approach for ten initial conditions of each model, ten spatially different stimuli, and different shock amplitudes. The investigated models of 2D excitable media cover a broad range of dominant frequencies and number of phase singularities, which demonstrates, that our findings are not limited to a specific kind of model or parameterization of it. Thus, we propose a method that incorporates the dynamics of the underlying system, even during pacing, and solely relies on a scalar observable, which is easily measurable in numerical simulations.

## 1 Introduction

Sudden cardiac death is the cause of 15%–20% of all deaths in western societies, making it a major health concern (Wong et al., 2019). Among the other types of cardiac arrhythmias, atrial fibrillation does not usually pose an immediate threat to life, but it is the most common form of sustained cardiac arrhythmia and often leads to strokes (Wolf et al., 1978). In contrast, ventricular fibrillation proves to be lethal if it is not stopped within minutes and is the most common cause of death in industrialised countries (Myerburg et al., 1997). Cardiac arrhythmias in network physiology are caused by chaotic patterns of the electrical excitation wave dynamics in the myocardium, that are not yet fully understood (Gray et al., 1998). Terminating cardiac arrhythmias like atrial and ventricular fibrillation is an ambitious task and most ofen requires the application of a high-energy electric shock (Ammannaya, 2020; Cheskes et al., 2022). As considerable electrical currents are induced in the myocardium, this approach involves severe side effects like post-traumatic stress (Godemann et al., 2004; Sears et al., 2011) and an increasing severity of life-threatening myocardial disorders after resuscitation (Xie et al., 1997; Tang et al., 2000).

We refer to defibrillation/pacing as the application of controlled electrical stimuli through an external electric field to counter spiral/scroll wave dynamics in the electrical excitation of the myocardium (Davidenko et al., 1995) with the ultimate goal to terminate all chaotic patterns in the heart rhythm. After successful termination of all chaotic patterns of the transmembrane potential of the heart, the sinus rhythm can be reinitiated by the electrical system of the heart.

There are numerous protocols for the control of chaotic patterns in excitable media with the aim of reducing the energy required for successful defibrillation. It is possible to apply stimuli at certain locations in the heart to terminate all chaotic patterns. For example, DeTal et al., 2022 identified specific isolines of the refractory backs of spiral waves as vulnerable locations. If stimuli are applied solely at these locations, the underlying dynamics is terminated. Lilienkamp and Parlitz, 2020 demonstrated that, in principle, it is sufficient to apply stimuli at isolated spots for successful termination. Otani et al., 2019 performed 3D simulations on a simplified ventricular geometry, to show the impact of the direction of an external electrical field on the success rate for defibrillation. Other protocols aim to apply either single stimuli at specific points in time (Trayanova, 2006; Steyer et al., 2023) or pulse sequences at predefined periods/frequencies (Fenton et al., 2009; Janardhan et al., 2014; Buran et al., 2017; Hornung et al., 2017; Lilienkamp et al., 2022a; Lilienkamp et al., 2022b; Aron et al., 2023). The timing of pulses can also be determined by feedback-controlled protocols, which incorporate the underlying dynamics and the reaction of the system to stimuli. Possible observables suitable for this group of algorithms were discussed by Buran et al., 2023 whereby the length of the refractory boundary has already been tested as a possible observable for a feedback-controlled pacing strategy (Buran, 2023). Numerical investigations of all these approaches promise reductions of the energy required for successful defibrillation and thus mitigate the aforementioned negative side effects of a conventional defibrillation. However, it remains unclear whether further energy reductions are in principle possible.

Numerical simulations of defibrillation attempts vary drastically in the level of detail of the underlying model. Striving to achieve a patient-specific level of parameterization, numerical simulations can be tailored accordingly, which involves the usage of the patients real heart geometry (Trayanova and Rantner, 2014; Lee et al., 2018), the orientation of its muscle fibres (Helm et al., 2005; Bayer et al., 2012), the cardiac cell model (ten Tusscher et al., 2004; Bueno-Orovio et al., 2008), etc. This type of modelling aims to predict patient-specific results (Niederer et al., 2019; Gerach et al., 2021; Gillette et al., 2021). However, multiple groups develop and test new pacing protocols on simplified models to express the temporal evolution of the transmembrane potential in 2D without the incorporation of the fibre orientation (Buran et al., 2017; Lilienkamp and Parlitz, 2020; DeTal et al., 2022). The low computation effort of running these simplified simulations allows for a flexible exploration of novel methods across a broad range of parameterizations, encompassing diverse models and parameterizations of them. Basic mechanisms for successful defibrillation can be identified but the accuracy of the predictions (e.g., the exact energy reduction) with respect to *in vivo* or *ex vivo* experiments is limited.

In this study, we investigate a novel feedback-based approach that incorporates the dynamics of the heart during pacing, to further minimize the energy required for successful defibrillation and thus mitigate the negative side effects associated with high-energy defibrillation. Compared to the aforementioned pacing protocols, the pulse timings derived with this approach are therefore not predetermined but calculated during the pacing process, based on the response of the system to the previous stimuli. The proposed approach solely relies on the measurement of the mean value of the transmembrane potential, an observable that can be measured very easily in numerical simulations.

In Section 2, we introduce the modelling and numerical treatment of our control attempts as well as a novel pacing protocol. A statistical analysis of the amplitude reduction achievable with our control scheme based on multiple cardiac tissue models, a comparison of its working mechanisms with pre-existing approaches, and a sensitivity analysis are given in Section 3.

## 2 Materials and methods

In this section we present the modelling and numerical treatment of our control attempts as well as a novel pacing protocol.

### 2.1 Simulation of cardiac dynamics

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

#### 2.1.1 Basic equations and numerical setup

The underlying mono-domain equations to model the temporal evolution of the transmembrane potential *V*_{m} are given by (Tuckwell, 1988)

where *D* describes the diffusion tensor, *I*_{ion} the ionic currents through the cell membrane, *C*_{m} the membrane capacitance per unit membrane area, and *I*_{ion}, *g*, and *h* depends on the model of cardiac dynamics in use, as described in Section 2.1.2 and the Supplementary Material. For all simulations we used no-flux boundary conditions

(where *n* describes the unit vector normal to the boundaries) and assumed isotropic as well as homogeneous diffusion *x* and *N* × *N* nodes using the finite difference method and fully explicit Euler method with Δ*t* (Press, 2007). To increase the numerical stability of the numerical scheme, we applied the Rush-Larsen (Rush and Larsen, 1978) method for the temporal integration of the gating variables. The chosen parameters Δ*x*, *N* × *N*, Δ*t*, and *D* meet the stability criterion of the FTCS scheme (Press (2007)) and were specified individually for each cardiac cell model to meet characteristic features (e.g., number of phase singularities) of the resulting chaotic states and can be found in Table 1.

**Table 1**. Modelling parameters for each cardiac cell model. Chaotic states can be characterized by the dominant frequency *f*_{dom} and the mean number of phase singularities *N*_{PS}. Electrical stimuli of length *t*_{pulse} are modelled via the concept of AVEs with a size *S*_{ave} of a single AVE, where a number *N*_{ave} of AVEs cover a relative fraction *A*_{cov} of the domain. Δ*t*_{median} is the median time between a local minimum and the consecutive local maximum in the mean value of the transmembrane potential.

#### 2.1.2 Investigated models of cardiac dynamics

Models of cardiac dynamics describe the electrical action potential dynamics of cardiomyocytes (Clayton et al., 2011). Different models of cardiac dynamics can be used to model the electrical activity in the myocardium of different species or different types of arrhythmia and thus exhibit different properties in the underlying dynamics, such as the shape and duration of action potentials (Fenton and Karma, 1998; Bueno-Orovio et al., 2008; Mahajan et al., 2008). We carried out simulations with four different models: the Aliev-Panfilov model (AP) (Aliev and Panfilov, 1996), the Bueno-Orovio-Cherry-Fenton model (BOCF) (Bueno-Orovio et al., 2008), the Fenton-Karma model (FK) (Fenton and Karma, 1998), and the Mitchell-Schaeffer model (Mitchell and Schaeffer, 2003) with the adaptions made by Alvarez et al. (MSA) (Álvarez et al., 2012). We specifically included the latter model because it represents an exception as it models ischaemic cardiac tissue. Although the other models can reconstruct distinctive features of diverse arrhythmia types, their parameterization is primarily designed to align with the shape and duration of the action potential during a state of sinus rhythm. The models we used for our experiments differ from more complex ionic models (ten Tusscher et al., 2004; Luo and Rudy, 1994) by not modelling the behaviour of all ion channel currents individually, but by combining the transmembrane currents/ion concentrations in specified functions/gating variables. Our choice of models and their parameterizations (given in the Supplementary Material) result in a broad range of characteristic features of the underlying dynamics, such as the number of phase singularities and dominant frequencies (see Table 1). This is intended to show that 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 mechanism inherent in all of these models. Due to a lower computational effort, the usage of phenomenological models allows us to perform a large number of simulations and thus a statistical analysis of the presented methods on multiple models of cardiac dynamics. Exemplary snapshots of the transmembrane potential *V*_{m} during spiral wave chaos are shown in Figure 1 for each model, respectively. As done in previous publications on pacing protocols, we also neglect the influence of the fibre orientation in the myocardium, assume absence of scar tissue and except for the MSA model, healthy electrophysiological behaviour of the ion channels within the myocardium (Trayanova, 2006; Janardhan et al., 2014; Buran et al., 2017; Hornung et al., 2017; Lilienkamp and Parlitz, 2020; Lilienkamp et al., 2022a; Lilienkamp et al., 2022b; DeTal et al., 2022; Steyer et al., 2023).

**Figure 1**. Snapshots of the membrane potential *V*_{m} during spiral wave chaos of **(A)** the Aliev-Panfilov (AP) model, **(B)** the Bueno-Orovio-Cherry-Fenton (BOCF) model, **(C)** the Fenton-Karma (FK) model, and **(D)** the Mitchell-Schaeffer model with the adaptions made by Alvarez et al. (MSA).

#### 2.1.3 Initialization of spiral wave chaos

For each model of cardiac dynamics, ten initial conditions were created by applying spatially randomized stimulations to the transmembrane potential of a state already exhibiting a spiral wave pattern. After these stimulations, each state was evolved for at least 5 s to generate independent initial conditions. Each initial condition was tested for another 5 s, to exclude self-terminating states (Qu, 2006; Lilienkamp et al., 2017; Aron et al., 2019). This way, we ensure that successful termination attempts can be attributed to the corresponding pacing protocol. Characteristics of the underlying dynamics created this way are given in Table 1. Further details about the initialization of chaotic states and video material of it are given in the Supplementary Material.

#### 2.1.4 Termination attempts

We simplify the application of electrical stimuli by using the concept of artificial virtual electrodes: 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). This underlying mechanism is modelled by the application of local current injections at predefined locations within the simulation domain (Lilienkamp et al., 2022b). These local current injections act as a network of local stimulation sites and are referred to as artificial virtual electrodes (AVEs). For each initial condition and each termination attempt, a fixed number *N*_{ave} of non-overlapping locations, of size *S*_{ave} are selected such that a relative fraction *A*_{cov} of the simulation domain is covered by local current injections. In previous publications, the coverage area *A*_{cov} was chosen from a broad range from 2.5% to 5% (Buran et al., 2017) to 25%–30% (Lilienkamp et al., 2022b) to 100% (Gray, 2002). The choice of *A*_{cov} is determined not only by the spatial resolution of the simulation domain but also by the decision which heterogeneities in the myocardium are modelled as virtual electrodes. We follow the abstraction of Lilienkamp et al., 2022b that different types of heterogeneities may contribute to the entirety of AVEs: blood vessels (Buran et al., 2017), lateral intercellular couping as well as extracellular resistive heterogeneities (Caldwell et al., 2015), boundaries of high curvature (Bittihn et al., 2012), and others.

The number of AVEs *N*_{ave}, the size of each AVE *S*_{ave}, and the resulting spatial fraction *A*_{cov} for each model of cardiac dynamics is given in Table 1. However, results for differenet values of *A*_{cov} are also presented in chapter 3.4, to emphasize the stability of our conclusion with respect to the exact choice of *A*_{cov}. An exemplary excitation using an AVE mask is shown in Figure 3A. The locally injected current of each AVE was selected randomly from the interval [0, *A*_{max}] to account for variations of depolarization intensity in living tissue due to the shape of heterogeneities. Relatively small heterogeneities, for example, are stimulated only at large external field strengths (Bittihn et al., 2012). By employing the concept of artificial virtual electrodes, we accept that in particular the first milliseconds of the interaction between an external electric field and the myocardium is simplified. 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. The described methodology matches the one used by Lilienkamp et al., 2022b. To determine the success rate of a pacing protocol at a given amplitude *A*_{max}, we performed ten termination attempts (each with different spatial distributions of AVEs) for each of the ten initial conditions. The success rate describes the fraction of successful attempts. A termination attempt is considered successful if after a period of 5˜*f*_{dom}, the maximum of the transmembrane potential is below a threshold *V*_{m,thresh} = 0.01 [*a*.*u*.]. The dominant frequency *f*_{dom} is defined as the frequency with the largest amplitude in the power spectral density of the underlying dynamics.

### 2.2 Local minima pacing

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}

as the observable to calculate pulse timings. Figure 2A illustrates, that local minima of *t*_{up} as the time between a local minimum and the consecutive local maximum in

**Figure 2**. Exemplary mean transmembrane potential time series of unperturbed (**(A)**, **(B)**) chaotic states for the FK model to illustrate the parameter definitions of Local Minima Pacing (LMP). All time series are scaled to [0,1] for ease of interpretation. **(A)** Time series of the mean transmembrane potential (grey) and mean values of the gating variables (yellow and cyan). Local minima in the mean value of the transmembrane potential **(B)** The upward deflection time Δ*t*_{up} (green) past local minima (yellow dots) in

In Local Minima Pacing (LMP), a pulse is applied shortly after a local minimum of *∀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 *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

**Algorithm 1.**Local Minima Pacing - Pseudocode.

**Input:**

**Input:** Δ*t*_{median} Median upward deflection time.

**Input:** *a* Scaling factor.

**Input:** *t*_{pause} Pause between two consecutive pulses.

**Input:** *max*_{pulses} maximum number of pulses to apply.

*num*_{p} ← 0, *t*_{last,p} ← 0

**while** not successful AND *num*_{p} < *max*_{pulses} **do**

**if** *t*_{c} − *t*_{last,p} > *t*_{pause} **then**

**if** *∀t* ∈ [*t*_{c} − *a* Δ*t*_{median}, *t*_{c}] **then**

Apply pulse

*t*_{last,p} ← *t*_{c}

*num*_{p} ← *num*_{p} + 1

**end if**

**end if**

*t*_{c} ← *t*_{c} + Δ*t*

**end while**

## 3 Results

In this section we present the outcomes derived with the Local Minima Pacing protocol, an investigation of optimal parameterisations of LMP for different models of cardiac dynamics, a comparison of the derived pulse timings with existing protocols, and a sensitivity analysis of the coverage rate *A*_{cov}.

### 3.1 Applying pulses past local minima

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

**Figure 3**. Exemplary termination attempt with Local Minima Pacing (LMP) for the Fenton-Karma model. **(A)** Snapshots of the transmembrane potential *V*_{m} at different points in time, corresponding to the purple vertical lines in subplot **(B)** where time series of the mean value of the transmembrane potential of the unperturbed (grey) and perturbed (green) dynamics are shown. LMP pulses (red lines) are applied shortly (0.2 Δ*t*_{median}) after local minima are detected in

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.

**Figure 4**. Dose-response curves for conventional defibrillation (grey squares), Equidistant Pacing (orange crosses), Adaptive Deceleration Pacing (blue triangles), and Local Minima Pacing (green circles) for **(A)** the Aliev-Panfilov (AP) model, **(B)** the Bueno-Orovio-Cherry-Fenton (BOCF) model, **(C)** the Fenton-Karma (FK) model, and **(D)** the Mitchell-Schaeffer model with the adaptions made by Alvarez et al. (MSA). The system was observed for 5 seconds to determine the pacing frequencies for Equidistant and Adaptive Deceleration Pacing. For Local Minima Pacing, the system was observed for only 2 seconds and the temporal distance between a local minimum and the application of a pulse was set to 0.2 Δ*t*_{median}.

### 3.2 Systematic investigation of pulse timings past local minima

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

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.

**Figure 5**. Systematic investigation of pulse timings for the Local Minima Pacing. Results for different pulse timings with respect to the required amplitude *A*_{max} to obtain a 90% success rate (colorbars) are shown for **(A)** the Aliev-Panfilov (AP) model, **(B)** the Bueno-Orovio-Cherry-Fenton (BOCF) model, **(C)** the Fenton-Karma (FK) model, and **(D)** the Mitchell-Schaeffer model with the adaptions made by Alvarez et al. (MSA), respectively. Grey squares indicate that *E*_{90} could not be achived with the given values of Δ*t*_{start} and Δ*t*_{end}. Purple crosses mark combinations of Δ*t*_{start} and Δ*t*_{end} for which the required E_{90}-energy of LMP is lower compared to all other protocols under investigation.

### 3.3 Comparison of LMP pulse timings with other protocols

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 *a* = 0.2, performed an *a posteriori* analysis of the temporal distances (periods) between LMP pulses at the respective E_{90}-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.

**Figure 6**. Analysis of pulse periods (temporal distance between two consecutive pulses) for the EquiDist, ADP, and LMP pacing protocols. Each column of a single plot denotes the pulse pair within a single pacing protocol, while the relative occurence of a given temporal distance (*y*-axis) over all initial conditions and all termination attempts is color coded. Plots **(A)**, **(D)**, **(G)**, and **(J)** reflect the distribution of 0.8 *f*_{dom}. Plots **(B)**, **(E)**, **(H)**, and **(K)** show the distribution of periods generated using the ADP approach. Plots **(C)**, **(F)**, **(I)**, and **(L)** visualize the *a posteriori* calculated LMP periods derived with the E90-amplitude. Occurrences of 0 ms at given pulse pairs indicate successful termination without applying further pulses.

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

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 *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.

### 3.4 Sensivity analysis of the modelled coverage area

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.

**Figure 7**. Effect of lower (marker: circle, triangle, cross) coverage rates *A*_{cov} on dose response curves of LMP (green) and ADP (blue) for **(A)** the Aliev-Panfilov (AP) model, **(B)** the Bueno-Orovio-Cherry-Fenton (BOCF) model, **(C)** the Fenton-Karma (FK) model, and **(D)** the Mitchell-Schaeffer model with the adaptions made by Alvarez et al.(MSA), respectively. Lower coverage areas require higher amplitudes *A*_{max} to achieve similar success rates.

## 4 Discussion

In this section, we outline the limitations and constraints of our study, highlight areas where further research is needed, and summarise our findings and their implications with respect to the heart as an essential part of the cardiovascular network.

### 4.1 Limitations

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.

### 4.2 Summary

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.

## Data availability statement

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

## Author contributions

DS: Conceptualization, Investigation, Software, Visualization, Writing–original draft, Writing–review and editing. SL: Writing–review and editing. TL: Conceptualization, Supervision, Writing–review and editing, Funding acquisition.

## Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—538874043.

## Acknowledgments

The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR project fhnz. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG)—440719683. The authors would like to thank the library of the Nuremberg Institute of Technology Georg Simon Ohm for covering the publication costs.

## 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.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

## 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.2024.1401661/full#supplementary-material

## Abbreviations

ADP, adaptive deceleration pacing; AP, Aliev-Panfilov model; AVE, artificial virtual electrode; BOCF, Bueno-Orovio-Cherry-Fenton model; EquiDist, equidistant pacing; FK, Fenton-Karma model; LMP, local minima pacing; MSA, Mitchell-Schaeffer model with adaptions by Alvarez; AF, atrial fibrillation; VF, ventricular fibrillation.

## References

Aliev, R. R., and Panfilov, A. V. (1996). A simple two-variable model of cardiac excitation. *Chaos Solit. Fractals* 7, 293–301. doi:10.1016/0960-0779(95)00089-5

Álvarez, D., Alonso-Atienza, F., Rojo-Álvarez, J. L., García-Alberola, A., and Moscoso, M. (2012). Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: a model study. *Math. Comput. Model.* 55, 1770–1781. doi:10.1016/j.mcm.2011.11.025

Ambrosi, D., Arioli, G., Nobile, F., and Quarteroni, A. (2011). Electromechanical coupling in cardiac dynamics: the active strain approach. *SIAM J. Appl. Math.* 71, 605–621. doi:10.1137/100788379

Ammannaya, G. K. K. (2020). Implantable cardioverter defibrillators–the past, present and future. *Archives Med. Science-Atherosclerotic Dis.* 5, 163–170. doi:10.5114/amsad.2020.97103

Aron, M., Lilienkamp, T., Luther, S., and Parlitz, U. (2023). Optimising low-energy defibrillation in 2D cardiac tissue with a genetic algorithm. *Front Netw Physiol.* 3, 1172454. doi:10.3389/fnetp.2023.1172454

Aron, M., Herzog, S., Parlitz, U., Luther, S., and Lilienkamp, T. (2019). Spontaneous termination of chaotic spiral wave dynamics in human cardiac ion channel models. *PloS one* 14, e0221401. doi:10.1371/journal.pone.0221401

Bayer, J. D., Blake, R. C., Plank, G., and Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. *Ann. Biomed. Eng.* 40, 2243–2254. doi:10.1007/s10439-012-0593-5

Bear, L. R., Dogrusoz, Y. S., Svehlíková, J., Coll-Font, J., Good, W., van Dam, E., et al. (2018). Effects of ecg signal processing on the inverse problem of electrocardiography. In *Comput. Cardiol., 2018 computing in cardiology conference (CinC)* (IEEE), 1–4. doi:10.22489/CinC.2018.070

Bittihn, P., Hörning, M., and Luther, S. (2012). Negative curvature boundaries as wave emitting sites for the control of biological excitable media. *Phys. Rev. Lett.* 109, 118106. doi:10.1103/PhysRevLett.109.118106

Bollmann, A., Kanuru, N., McTeague, K., Walter, P., DeLurgio, D., and Langberg, J. (1998). Frequency analysis of human atrial fibrillation using the surface electrocardiogram and its response to ibutilide. *Am. J. Cardiol.* 81, 1439–1445. doi:10.1016/s0002-9149(98)00210-0

B Traykov, V., Pap, R., and Sághy, L. (2012). Frequency domain mapping of atrial fibrillation-methodology, experimental data and clinical implications. *Curr. Cardiol. Rev.* 8, 231–238. doi:10.2174/157340312803217229

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

Buran, C. P. (2023) *The mechanism of fibrillation control by low-energy anti-fibrillation pacing*, 70–107.

Buran, P., Bär, M., Alonso, S., and Niedermayer, T. (2017). Control of electrical turbulence by periodic excitation of cardiac tissue. *Chaos An Interdiscip. J. Nonlinear Sci.* 27, 113110. doi:10.1063/1.5010787

Buran, P., Niedermayer, T., and Baer, M. (2023). Mechanism of defibrillation of cardiac tissue by periodic low-energy pacing. *bioRxiv*, 2023–2103. doi:10.1101/2023.03.16.533010

Caldwell, B. J., Trew, M. L., and Pertsov, A. M. (2015). Cardiac response to low-energy field pacing challenges the standard theory of defibrillation. *Circulation Arrhythmia Electrophysiol.* 8, 685–693. doi:10.1161/CIRCEP.114.002661

Chen, J., Mandapati, R., Berenfeld, O., Skanes, A. C., and Jalife, J. (2000). High-frequency periodic sources underlie ventricular fibrillation in the isolated rabbit heart. *Circulation Res.* 86, 86–93. doi:10.1161/01.res.86.1.86

Cheskes, S., Verbeek, P. R., Drennan, I. R., McLeod, S. L., Turner, L., Pinto, R., et al. (2022). Defibrillation strategies for refractory ventricular fibrillation. *N. Engl. J. Med.* 387, 1947–1956. doi:10.1056/NEJMoa2207304

Clayton, R., Bernus, O., Cherry, E., Dierckx, H., Fenton, F. H., Mirabella, L., et al. (2011). Models of cardiac tissue electrophysiology: progress, challenges and open questions. *Prog. biophysics Mol. Biol.* 104, 22–48. doi:10.1016/j.pbiomolbio.2010.05.008

Davidenko, J. M., Salomonsz, R., Pertsov, A. M., Baxter, W. T., and Jalife, J. (1995). Effects of pacing on stationary reentrant activity: theoretical and experimental study. *Circulation Res.* 77, 1166–1179. doi:10.1161/01.res.77.6.1166

DeTal, N., Kaboudian, A., and Fenton, F. H. (2022) “Terminating spiral waves with a single designed stimulus: teleportation as the mechanism for defibrillation,”in, *Proceedings of the national academy of sciences*. doi:10.1073/pnas.2117568119

Efimov, I. R., Aguel, F., Cheng, Y., Wollenzier, B., and Trayanova, N. (2000). Virtual electrode polarization in the far field: implications for external defibrillation. *Am. J. Physiology-Heart Circulatory Physiology* 279, H1055–H1070. doi:10.1152/ajpheart.2000.279.3.H1055

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

Fenton, F. H., Luther, S., Cherry, E. M., Otani, N. F., Krinsky, V., Pumir, A., et al. (2009). Termination of atrial fibrillation using pulsed low-energy far-field stimulation. *Circulation* 120, 467–476. doi:10.1161/CIRCULATIONAHA.108.825091

Gerach, T., Schuler, S., Fröhlich, J., Lindner, L., Kovacheva, E., Moss, R., et al. (2021). Electro-mechanical whole-heart digital twins: a fully coupled multi-physics approach. *Mathematics* 9, 1247. doi:10.3390/math9111247

Gillette, K., Gsell, M. A., Prassl, A. J., Karabelas, E., Reiter, U., Reiter, G., et al. (2021). A framework for the generation of digital twins of cardiac electrophysiology from clinical 12-leads ecgs. *Med. Image Anal.* 71, 102080. doi:10.1016/j.media.2021.102080

Godemann, F., Butter, C., Lampe, F., Linden, M., Schlegl, M., Schultheiss, H.-P., et al. (2004). Panic disorders and agoraphobia: side effects of treatment with an implantable cardioverter/defibrillator. *Clin. Cardiol.* 27, 321–326. doi:10.1002/clc.4960270604

Gray, R. A. (2002). Termination of spiral wave breakup in a fitzhugh–nagumo model via short and long duration stimuli. *Chaos An Interdiscip. J. Nonlinear Sci.* 12, 941–951. doi:10.1063/1.1497836

Gray, R. A., and Chattipakorn, N. (2005). Termination of spiral waves during cardiac fibrillation via shock-induced phase resetting. *Proc. Natl. Acad. Sci. U. S. A.* 102, 4672–4677. doi:10.1073/pnas.0407860102

Gray, R. A., Pertsov, A. M., and Jalife, J. (1998). Spatial and temporal organization during cardiac fibrillation. *Nature* 392, 75–78. doi:10.1038/32164

Hazim, A., Belhamadia, Y., and Dubljevic, S. (2021). A simulation study of the role of mechanical stretch in arrhythmogenesis during cardiac alternans. *Biophysical J.* 120, 109–121. doi:10.1016/j.bpj.2020.11.018

Helm, P., Beg, M. F., Miller, M. I., and Winslow, R. L. (2005). Measuring and mapping cardiac fiber and laminar architecture using diffusion tensor mr imaging. *Ann. N. Y. Acad. Sci.* 1047, 296–307. doi:10.1196/annals.1341.026

Hooks, D. A., Trew, M. L., Caldwell, B. J., Sands, G. B., LeGrice, I. J., and Smaill, B. H. (2007). Laminar arrangement of ventricular myocytes influences electrical behavior of the heart. *Circulation Res.* 101, e103–e112. doi:10.1161/CIRCRESAHA.107.161075

Hornung, D., Biktashev, V., Otani, N., Shajahan, T., Baig, T., Berg, S., et al. (2017). Mechanisms of vortices termination in the cardiac muscle. *R. Soc. open Sci.* 4, 170024. doi:10.1098/rsos.170024

Hu, Y., Gurev, V., Constantino, J., Bayer, J. D., and Trayanova, N. A. (2013). Effects of mechano-electric feedback on scroll wave stability in human ventricular fibrillation. *PloS one* 8, e60287. doi:10.1371/journal.pone.0060287

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

Leclercq, C., Faris, O., Tunin, R., Johnson, J., Kato, R., Evans, F., et al. (2002). Systolic improvement and mechanical resynchronization does not require electrical synchrony in the dilated failing heart with left bundle-branch block. *Circulation* 106, 1760–1763. doi:10.1161/01.cir.0000035037.11968.5c

Lee, A. W., Costa, C. M., Strocchi, M., Rinaldi, C. A., and Niederer, S. A. (2018). Computational modeling for cardiac resynchronization therapy. *J. Cardiovasc. Transl. Res.* 11, 92–108. doi:10.1007/s12265-017-9779-4

Lilienkamp, T., Christoph, J., and Parlitz, U. (2017). Features of chaotic transients in excitable media governed by spiral and scroll waves. *Phys. Rev. Lett.* 119, 054101. doi:10.1103/PhysRevLett.119.054101

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

Lilienkamp, T., Parlitz, U., and Luther, S. (2022a). Non-monotonous dose response function of the termination of spiral wave chaos. *Sci. Rep.* 12, 12043. doi:10.1038/s41598-022-16068-8

Lilienkamp, T., Parlitz, U., and Luther, S. (2022b). Taming cardiac arrhythmias: terminating spiral wave chaos by adaptive deceleration pacing. *Chaos An Interdiscip. J. Nonlinear Sci.* 32, 121105. doi:10.1063/5.0126682

Luo, C.-h., and Rudy, Y. (1994). A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes. *Circulation Res.* 74, 1071–1096. doi:10.1161/01.res.74.6.1071

Luther, S., Fenton, F. H., Kornreich, B. G., 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

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L.-H., et al. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. *Biophysical J.* 94, 392–410. doi:10.1529/biophysj.106.98160

Mitchell, C. C., and Schaeffer, D. G. (2003). A two-current model for the dynamics of cardiac membrane. *Bull. Math. Biol.* 65, 767–793. doi:10.1016/S0092-8240(03)00041-7

Myerburg, R. J., Interian Jr, A., Mitrani, R. M., Kessler, K. M., and Castellanos, A. (1997). Frequency of sudden cardiac death and profiles of risk. *Am. J. Cardiol.* 80, 10F-19F–19F. doi:10.1016/s0002-9149(97)00477-3

Nash, M. P., Mourad, A., Clayton, R. H., Sutton, P. M., Bradley, C. P., Hayward, M., et al. (2006). Evidence for multiple mechanisms in human ventricular fibrillation. *Circulation* 114, 536–542. doi:10.1161/CIRCULATIONAHA.105.602870

Niederer, S. A., Lumens, J., and Trayanova, N. A. (2019). Computational models in cardiology. *Nat. Rev. Cardiol.* 16, 100–111. doi:10.1038/s41569-018-0104-y

Otani, N. F., Wheeler, K., Krinsky, V., and Luther, S. (2019). Termination of scroll waves by surface impacts. *Phys. Rev. Lett.* 123, 068102. doi:10.1103/PhysRevLett.123.068102

Pathmanathan, P., Cordeiro, J. M., and Gray, R. A. (2019). Comprehensive uncertainty quantification and sensitivity analysis for cardiac action potential models. *Front. physiology* 10, 721. doi:10.3389/fphys.2019.00721

Press, W. H. (2007). “Numerical recipes,” in *The art of scientific computing*. 3rd edition (New York, NY, United States: Cambridge University Press).

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

Pumir, A., Nikolski, V., Hörning, M., Isomura, A., Agladze, K., Yoshikawa, K., et al. (2007). Wave emission from heterogeneities opens a way to controlling chaos in the heart. *Phys. Rev. Lett.* 99, 208101. doi:10.1103/PhysRevLett.99.208101

Qu, Z. (2006). Critical mass hypothesis revisited: role of dynamical wave stability in spontaneous termination of cardiac fibrillation. *Am. J. Physiology-Heart Circulatory Physiology* 290, H255–H263. doi:10.1152/ajpheart.00668.2005

Ripplinger, C. M., Krinsky, V. I., Nikolski, V. P., and Efimov, I. R. (2006). Mechanisms of unpinning and termination of ventricular tachycardia. *Am. J. Physiology-Heart Circulatory Physiology* 291, H184–H192. doi:10.1152/ajpheart.01300.2005

Rush, S., and Larsen, H. (1978). A practical algorithm for solving dynamic membrane equations. *IEEE Trans. Biomed. Eng.* 25, 389–392. doi:10.1109/TBME.1978.326270

Salvador, M., Regazzoni, F., Dede, L., Quarteroni, A., et al. (2023). Fast and robust parameter estimation with uncertainty quantification for the cardiac function. *Comput. Methods Programs Biomed.* 231, 107402. doi:10.1016/j.cmpb.2023.107402

Sears, S. F., Hauf, J. D., Kirian, K., Hazelton, G., and Conti, J. B. (2011). Posttraumatic stress and the implantable cardioverter-defibrillator patient: what the electrophysiologist needs to know. *Circulation Arrhythmia Electrophysiol.* 4, 242–250. doi:10.1161/CIRCEP.110.957670

Shajahan, T., Berg, S., Luther, S., Krinski, V., and Bittihn, P. (2016). Scanning and resetting the phase of a pinned spiral wave using periodic far field pulses. *New J. Phys.* 18, 043012. doi:10.1088/1367-2630/18/4/043012

Skanes, A. C., Mandapati, R., Berenfeld, O., Davidenko, J. M., and Jalife, J. (1998). Spatiotemporal periodicity during atrial fibrillation in the isolated sheep heart. *Circulation* 98, 1236–1248. doi:10.1161/01.cir.98.12.1236

Stewart, A., Allen, J., and Adgey, A. (1992). Frequency analysis of ventricular fibrillation and resuscitation success. *QJM An Int. J. Med.* 85, 761–769. doi:10.1093/oxfordjournals.qjmed.a068713

Steyer, J., Lilienkamp, T., Luther, S., and Parlitz, U. (2023). The role of pulse timing in cardiac defibrillation. *Front. Netw. Physiology* 2, 1007585. doi:10.3389/fnetp.2022.1007585

Tang, W., Weil, M. H., and Sun, S. (2000). Low-energy biphasic waveform defibrillation reduces the severity of postresuscitation myocardial dysfunction. *Crit. Care Med.* 28, N222–N224. doi:10.1097/00003246-200011001-00014

Ten Tusscher, K. H., Hren, R., and Panfilov, A. V. (2007). Organization of ventricular fibrillation in the human heart. *Circulation Res.* 100, e87–e101. doi:10.1161/CIRCRESAHA.107.150730

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

Trayanova, N. (2006). Defibrillation of the heart: insights into mechanisms from modelling studies. *Exp. Physiol.* 91, 323–337. doi:10.1113/expphysiol.2005.030973

Trayanova, N. A., and Rantner, L. J. (2014). New insights into defibrillation of the heart from realistic simulation studies. *Europace* 16, 705–713. doi:10.1093/europace/eut330

Tuckwell, H. C. (1988) *Cambridge studies in mathematical biology introduction to theoretical neurobiology: series number 8: linear cable theory and dendritic structure volume 1*. Cambridge, England: Cambridge University Press.

Wolf, P. A., Dawber, T. R., Thomas Jr, H. E., and Kannel, W. B. (1978). Epidemiologic assessment of chronic atrial fibrillation and risk of stroke: the fiamingham study. *Neurology* 28, 973. doi:10.1212/wnl.28.10.973

Wong, C. X., Brown, A., Lau, D. H., Chugh, S. S., Albert, C. M., Kalman, J. M., et al. (2019). Epidemiology of sudden cardiac death: global and regional perspectives. *Heart, Lung Circulation* 28, 6–14. doi:10.1016/j.hlc.2018.08.026

Xie, J., Weil, M. H., Sun, S., Tang, W., Sato, Y., Jin, X., et al. (1997). High-energy defibrillation increases the severity of postresuscitation myocardial dysfunction. *Circulation* 96, 683–688. doi:10.1161/01.cir.96.2.683

Keywords: defibrillation, arrhythmias, chaos control, excitable media, numerical simulation, feedback-controlled protocol, network physiology, cardiovascular network

Citation: Suth D, Luther S and Lilienkamp T (2024) Chaos control in cardiac dynamics: terminating chaotic states with local minima pacing. *Front. Netw. Physiol.* 4:1401661. doi: 10.3389/fnetp.2024.1401661

Received: 15 March 2024; Accepted: 26 April 2024;

Published: 03 July 2024.

Edited by:

Alessio Gizzi, Campus Bio-Medico University, ItalyReviewed by:

Simone Scacchi, University of Milan, ItalyRichard H. Clayton, The University of Sheffield, United Kingdom

Richard A. Gray, United States Food and Drug Administration, United States

Copyright © 2024 Suth, Luther and Lilienkamp. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Daniel Suth, daniel.suth@th-nuernberg.de