Optimising low-energy defibrillation in 2D cardiac tissue with a genetic algorithm

Sequences of low-energy electrical pulses can effectively terminate ventricular fibrillation (VF) and avoid the side effects of conventional high-energy electrical defibrillation shocks, including tissue damage, traumatic pain, and worsening of prognosis. However, the systematic optimisation of sequences of low-energy pulses remains a major challenge. Using 2D simulations of homogeneous cardiac tissue and a genetic algorithm, we demonstrate the optimisation of sequences with non-uniform pulse energies and time intervals between consecutive pulses for efficient VF termination. We further identify model-dependent reductions of total pacing energy ranging from ∼4% to ∼80% compared to reference adaptive-deceleration pacing (ADP) protocols of equal success rate (100%).

To address this need, pacing methods have been developed that aim to replace the single high-energy shock with a sequence of low-energy pulses. In general, these control algorithms differ in detail in the generation of the pulse sequences and are either open or closed-loop. Low-energy anti-fibrillation pacing (LEAP) uses a sequence of electrical pulses of constant amplitude, duration, and period to terminate atrial and ventricular fibrillation (Fenton et al., 2009;Luther et al., 2011;Buran et al., 2022). Using this approach in in-vivo experiments, a pulse energy reduction of 80%-90% compared to the conventional single shock has been achieved. However, several studies showed that the performance of pacing protocols with equidistant pulses depends very sensitively on the choice of the pacing interval (Buran et al., 2017;DeTal and Fenton, 2021;Lilienkamp et al., 2022a;Buran et al., 2023). As another example, multi-stage pacing developed by the Efimov lab employs groups of pulses with constant amplitude and timing (Janardhan et al., 2014).
Experimentally, optogenetically modified cardiac cells offer new ways for implementing and testing alternative approaches for defibrillation (Bruegmann et al., 2016;Hussaini et al., 2021). Progress has also been been made combining experimental data with numerical modelling, for example, for studying spatiotemporal alternans patterns of cardiac excitation (Loppini et al., 2022). Recent numerical studies have shown that pulse timing is crucial for termination success (Steyer et al., 2023). Furthermore, chaotic cardiac arrhythmias may exhibit self-termination (Rappel et al., 2022), a phenomenon also seen with chaotic transients in generic excitable media (Lilienkamp et al., 2017;Lilienkamp and Parlitz, 2018;Aron et al., 2019). The fact that transient chaos in excitable media can in principle be terminated by very few local perturbations (Lilienkamp and Parlitz, 2020) provides additional evidence that low-energy defibrillation is in principle possible. And in fact, recent studies show progress in the quest of finding better defibrillation methods, for example, exploiting targeted manipulation of spiral waves (DeTal et al., 2022), non-monotonous dose-response curves (Lilienkamp et al., 2022a) and decelerated pulse trains (Lilienkamp et al., 2022b).
In particular, the progress made with non-equidistant pulse sequences (Lilienkamp et al., 2022b) raises the question of further increasing termination efficiency when more complex sequences of (local) perturbations are used. However, the systematic experimental optimisation of such arrhythmia control algorithms with a large parameter space is under practical and ethical constraints not feasible. Therefore, in contrast to the above mechanistic, hypothesis-driven studies, we formulate the control of spatio-temporal dynamics in this study as an optimisation problem and use a genetic algorithm (GA) to cope with the vast landscape of possible pacing protocols spanning all combinations of pulse counts, amplitudes, and intervals between pulses.
Following natural evolution, the genetic algorithm mutates and combines pacing sequences and evaluates their "fitness" in terms of (high) termination probability and (low) pulse amplitudes (or energies). Because the goal of this study is to demonstrate the feasibility of GA optimisation of low-energy defibrillation in a proof of principle, we use simple phenomenological models of cardiac tissue in a homogeneous, isotropic 2D medium to avoid excessive computational burden that would result from more realistic modeling. This requires an appropriate selection of cardiac-tissue models as a foundation for simulations of twodimensional sheets of homogeneous cardiac tissue, as well as a simplified translation of real-life defibrillation into something a computer can parse in reasonable amounts of time. This simulation framework will be introduced in the following section.

Simulating ventricular fibrillation
In this section, we present the cardiac-tissue model(s) and numerical integration required to study and evaluate different pacing sequences.

Cardiac-tissue models
Cardiac-tissue models describe the electrophysiological dynamics of the heart. The heart contracts (i.e., pumps) in response to electrical excitation-waves propagating through its muscle tissue in a coordinated fashion. Such waves originate in the sinoatrial node in regular intervals to drive the heartbeat. Because cardiac tissue constitutes an excitable medium, orderly propagation is generally ensured: Each cell only responds to (and propagates) incoming signals above a certain activation threshold and then enters a temporary post-excitation refractory period where further propagation is greatly impeded. While the former prevents (weak) abnormal activation, the latter ensures the unidirectional movement of the signal for the coordinated contraction of each heart chamber. Both mechanisms are a result of ion channels permeating cardiac muscle cells, which open and close in a predetermined fashion upon activation. They allow for the migration of various ion species to and from the cell's interior, changing the net voltage across the cell membrane in the process before eventually returning it to its resting state.
We limited ourselves to phenomenological cardiac-tissue models for the pacing simulations for reasons of efficiency. Such models attempt to establish an abstract, high-level description of the complex physics involved in cardiac electrophysiology at manageable levels of mathematical complexity. This means they naturally lend themselves to large-scale computer simulations and experiments as they are usually both simpler to implement and less demanding of computational resources. These are in contrast to physiological models, which incorporate some (or more) of the underlying physiological details contributing to the signalpropagation dynamics in cardiac tissue. The incorporation of such details usually incurs significant additional computational cost, and we thus restricted our study to phenomenological models to reduce overall simulation time.
A generic cardiac-tissue model is given by a reaction-diffusion equation for the transmembrane potential u, the net voltage across a local aggregate (continuum) of cell membranes: A local net current I tot makes up the reaction term which may amplify or counteract the ion-diffusion process to neighbouring cells as governed by the diffusion tensor D. The net current consists of a model-dependent number of contributing currents modulated by a number of ion-channel gates whose permeabilities are represented by a set of dimensionless gating variables. C m is the cell-membrane capacitance. Lastly, the u is subject to a no-flux Neumann boundary condition where ∇u ·n 0 (withn as the boundary normal) holds everywhere on the physical boundary of the simulated tissue.
Throughout this paper, we generally refer to chaotic states of cardiac-tissue models as fibrillation episodes in reference to the actual medical condition of (ventricular) fibrillation which also features chaotic electrophysiological dynamics. The simulated control-methods used, which aim to terminate such chaotic states, are analogously called defibrillation methods in this context. However, it should be understood that two-dimensional simulations of homogeneous cardiac tissue and their control by local current injection imitating virtual electrodes represent a substantial simplification of the real heart. Accordingly, great care must be taken in generalising any observations and findings to the broader medical context.

The Fenton-Karma model
Our first model of choice was the Fenton-Karma (FK) model (Fenton and Karma (1998); Figure 1) for human ventricular tissue. It offers great versatility in its behaviour through its various high-level parameters (Fenton et al., 2002) and it can in fact be fine-tuned to approximate the dynamics of more complex (e.g., physiological) models as needed. Employing multiple FK parameter sets (see Tables 1, 2) enabled us to test any potential pacing-protocol improvements for one set for reproducibility with others at little extra computational cost.
The model itself features three dimensionless state variables: the transmembrane potential (u) and two gating variables (v, w). The gating variables govern two of the three transmembrane currents found in this model's net-current equation: where H denotes the standard Heaviside step-function and I stim is an external stimulation (e.g., pacing-induced) current.
is dependent on u, making it a time constant in need of constant re-evaluation during the integration process. The diffusion tensor D was assumed isotropic with a scalar value of D = 0.20 mm 2 /ms and the capacitance set to C m = 1 ms. The FK model's transmembrane voltage u is normalised such that its resting value lies at u = 0.0 a.u. for all its parameter sets, making comparisons between them straightforward.

The Bueno-Orovio-Cherry-Fenton model
The Bueno-Orovio-Cherry-Fenton (BOCF) model (Bueno-Orovio et al. (2008); Figure 2) was our second cardiac-tissue model of choice, which also describes the dynamics in human ventricular tissue. Like the FK model, the BOCF model is also capable of approximating other models on a qualitative level given proper parameters. It features four currents and three gating variables: The model has seven time constants, mostly in its τ parameters. Its resting potential is also normalised to 0.0 a.u. and the membrane capacitance set to unity (ms). We similarly assumed an isotropic D = 0.30 mm 2 /ms as well. We used a single parameter set (see Tables 2, 3) designed to emulate the behaviour of a much more complex physiological model in our simulations.

Numerical implementation
We limited simulations to square, isotropic 2D tissue sheets for speed and ease of implementation. We chose domain sizes spanning 256 2 and 512 2 pixels for the FK and BOCF model(s), respectively. Both had the same discretisation parameters of Δt = 0.1 ms for the time step and h def .

FIGURE 2
Snapshot of BOCF voltage dynamics. Area-wise, this domain is four times larger than the one used for the FK models, and the peak voltage of the model's action potential is also notably higher at approx. 1.5 a.u. in comparison.  We discretised voltage variables with a Forward in Time, Central in Space (FTCS) scheme, which is an explicit first-order method in both time and space (Li et al., 2018). For the Laplacian, we used a 9stencil approximation (Lynch, 1992) in place of the 5-stencil at nonvertex voltage-grid locations (i, j) to ensure a more radially independent discretisation error for proper spiral-wave propagation.
The separate treatment of vertices in the computation of the Laplacian is owed to our reliance on the "ghost point" method for the enforcement of the no-flux boundary condition on the voltage grid u (Li et al., 2018). This method extends the grid representing u by 2 pixels along each axis (one per side for a rectangular domain), creating a one-pixel boundary around it. Each such pixel is linked to another specific pixel of the interior grid (either one row or column apart) such that they both always hold the same value. Applying FTCS and Laplace discretisations on the main part of this extended grid then enforces the boundary condition automatically. The 5-stencil approximation is thus needed to compute the Laplacian on any of the interior grid's four vertices to avoid referencing the undefined values in any of the extension's four corners.
3 Modelling and evaluating pacing protocols 2D simulations of heterogeneous cardiac tissue by themselves lack the means to reproduce the defibrillation mechanisms of real hearts in an emergent fashion. We therefore had to model and introduce the qualitative properties of these mechanisms to our simulations manually to enable the simulation of pacing protocols. We focused on the mode of operation of low-energy pacing in particular for the modelling process.

Low-energy pacing and virtual electrodes
Despite advances in defibrillator design, the currents and energies involved still carry substantial risk to the patient's longterm health and may cause internal tissue damage or even long-term changes in a given patient's heart rhythm. This is the principal motivation behind contemporary low-energy defibrillation research (Pumir et al., 2007), which seeks to discover and study novel, potentially less intrusive pacing approaches such as low-energy anti-fibrillation pacing (LEAP) (Luther et al., 2011). LEAP replaces the singular, potent mono-or biphasic shocks with a sequence of temporally equidistant, uniformly low-amplitude ones to lower the risk of harmful side effects to the patient.
Low-energy pacing cannot rely on overwhelming the fibrillating signal in the heart by forcefully activating most muscle tissue; it exploits the numerous small conductivity heterogeneities (e.g., blood vessels, scar or fatty tissue) permeating cardiac tissue to perturb the electric activity in the heart instead. Such a heterogeneity can function as a local signal wave-front emitter by acting as an anode-cathode pair (a virtual electrode) when exposed to an electric field of sufficient strength and proper orientation. Such locally emitted wave fronts can then interact with the fibrillating signal in the heart by blocking its propagation through the refactory period they induce at various locations in the affected tissue, ideally putting an end to the fibrillation episode.
A virtual electrode's activation threshold and response depend on its geometry and orientation relative to the external electric field (Bittihn et al., 2008;Bittihn et al., 2012), where higher field strengths can generally be used to "recruit" more of these electrodes for defibrillation purposes. Low-energy defibrillation methods capitalise on low-threshold virtual electrodes which can trigger in response to weaker fields, potentially disrupting a fibrillation episode through a few periodic (in place of just a single potent one as used in traditional defibrillation) activations and subsequent emissions at minimal overall risk to the patient. While low-energy defibrillation methods such as LEAP have seen success in laboratory settings, they are still pending adoption in the medical context until further clinical trials have fully established their overall viability.

Our simplified pacing model
A pacing protocol comprises a sequence of pacing amplitudes and inter-pulse intervals. We preemptively reduced the dimensionality of the protocol-optimisation space by limiting the search to protocols of a uniform pulse length of 2.0 ms throughout all our simulations. Similarly, we further only considered monophasic protocols of positive amplitude(s) for our study. Under these constraints, a monophasic N-pulse protocol therefore featured N amplitudes and N − 1 inter-pulse intervals for a total of 2 N − 1 parameters. We conducted our investigation with 5-pulse protocols first, which still left us with 9 protocol parameters to optimise.
In a defibrillation attempt based on a given protocol, we would first apply it through randomly distributed injection sites modelling the "recruitment" of local virtual electrodes for pacing purposes as utilised in low-energy defibrillation approaches. We then observed the system for a short amount of time and employed a simple binary criterion based on the resulting peak transmembrane voltage across the tissue domain to check for successful defibrillation afterwards: If max(u(t)) < u crit. = 5 × 10 −2 a.u. holds across the entire tissue domain at time t, then no further action potentials can occur and the system is rapidly approaching its resting state. We used a global value for u crit. as all our cardiac-tissue models of choice featured normalised voltages variables.
We applied pacing protocols numerically through local, timedependent changes in a given cardiac-tissue model's stimulation current I stim ( x, t) as specified by protocol parameters. The individual areas of stimulation were chosen by randomly distributing non-overlapping injection sites spanning 3 × 3 pixels each (Figure 3) throughout the tissue domain. Each injection site was further associated with a uniformly drawn modulation coefficient between 0 and 1 which defined the (relative) local pacing efficacy at a given site to account for the effects of nonuniform virtual-electrode geometries and orientations in live tissue. We chose an injection-site coverage of~28% as an arbitrary starting point to account for the number of virtual electrodes involved, which translated to 2,000 and 8,000 injection sites for the FK and BOCF model(s), respectively. Once generated, both the injection-site locations and their respective modulation coefficients remained constant throughout the application of a given protocol.

Frontiers in Network Physiology
frontiersin.org The length of the post-pacing observation period was modeldependent at a value of 5 × T dom , where T dom denotes the dominant period of oscillation for a given model. While potentially chaotic systems like the cardiac-tissue models may not conform to a single period of oscillation, it is still possible to estimate a dominant (i.e., most prevalent in terms of its absolute Fourier-spectrum contribution) one by means of fast Fourier transforms (Schwaderlapp et al., 2017). This period can be understood as something resembling a mean period of oscillation for the wave fronts in the chaotic dynamics. As we generally want defibrillation to occur in a reasonable time frame, we settled on this 5 × T dom -interval so as to enforce this criterion.
To estimate the models' dominant periods (see Table 2), we integrated and observed one initial condition per model over an interval of 10 s at 120 samples per second, yielding one time-series of 1200 data points per pixel of the discretised voltage domain. The separation of the domain into individual time-series in this step can be justified with the spatiotemporally chaotic nature of the dynamics, where spatial (and temporal) correlations decay over short distances (or time intervals). We then applied a Fast Fourier Transform to the resulting time series, computed their highest-contributing frequencies, and subsequently averaged them to estimate "the" dominant frequency (and therefore period) of the given model.
We detected successful defibrillation following the post-pacing observation period by checking whether max(u) < 0.05 a.u. held for the resulting voltage grid. If it did, the system would be incapable of creating further wave fronts on its own and could therefore be considered defibrillated. We determined this threshold manually and applied it uniformly to all our selected cardiac-tissue models because of their common voltage ranges.

Defining protocol-performance metrics
We defined two performance metrics for pacing protocols: pacing cost (PC) and termination ratio (TR). PC measures the theoretical cost in applying a protocol based on its pulse amplitudes (analogous to energy; Figure 4), with lower values being preferred. A given protocol's TR is an estimate of its efficacy in terminating (i.e., defibrillating) simulated episodes of Sketch of the pacing model. Non-overlapping injection sites (green) are distributed across the voltage domain (right panel). Each site has a random efficacy coefficient to modulate the amplitude(s) of the pacing protocol to be applied (left panel). This emulates the effects of different heterogeneity geometries on a local level.

FIGURE 4
Sketch of the pacing-cost metric. The square of each pulse amplitude is integrated over its respective duration, and the sum over all pulses then yields the total cost. The squaring is done to penalise protocols with overly potent singular pulses.
Frontiers in Network Physiology frontiersin.org fibrillation; it ranges in value from 0 (completely ineffective) to 1 (peak efficacy). While PC can be computed at negligible cost for a given protocol, its TR value needs to be estimated by sampling the protocol's defibrillation efficacy through computationally expensive simulations. We conducted the sampling process through a series of defibrillation simulations involving a predetermined set of 10 model-dependent initial conditions (ICs) representing fibrillation states. This process was repeated a predetermined number of times to account for the influence of different injectionsite distributions (e.g., 5 distributions each, amounting to 5 × 10 samples per protocol). The TR value was then set to the resulting success (i.e., defibrillation) ratio among these samples. Two videos demonstrating a successful and failed termination attempt, respectively, are provided as Supplementary Material.
It should be noted that fibrillation episodes in our cardiac-tissue models of choice are of an inherently transient nature (Strain and Greenside, 1998;Lilienkamp et al., 2017;Lilienkamp and Parlitz, 2018;Aron et al., 2019) such that any chaotic trajectory ultimately converges to a steady (i.e., "defibrillated") state of its own, and this also applies to the general post-pacing (i.e., perturbation) state of such an episode. As the mean transient time for this transition grows exponentially with the physical domain size to which the dynamics are confined, our choice of sufficiently short post-pacing observation periods kept potential biases in our TR estimates through such "selfterminated" fibrillation episodes to a minimum.

Choosing reference protocols
We needed reference points to gauge prospective pacingoptimisation results. To that end, we consulted the adaptive (spectrum-guided) deceleration pacing (ADP) scheme (Lilienkamp et al., 2022b) to derive an "optimal" set of reference protocols.
ADP protocols feature pulses of uniform strength and duration in a decelerating fashion: the resting period between two consecutive excitations grows after each pulse. Such protocols serve as a rather strict point of comparison as they have been shown to outperform LEAP analogues at lower overall PC in simulations of homogeneous 2D cardiac tissue (Lilienkamp et al., 2022b).
The pulse timings are derived through the Fourier spectrum of a given cardiac-tissue model: a sufficiently long (10 s) timeseries of appropriate temporal resolution (10 ms) is first used to establish the spectrum and the cumulative integral of its squared absolute-value up to a given cutoff frequency. Partitioning this (normalised) integral equally according to the desired number of pulses then yields the corresponding inter-pulse periods ( Figure 5).
The resulting ADP protocol retains one degree of freedom in its uniform pulse-amplitudes. We thus evaluated the TR of each model's ADP protocol over multiple amplitude values to establish the corresponding PC thresholds for 90%, 95%, 99%, and 100% TR in a dose-response curve (Figure 6), yielding four reference protocols per model (Table 4).

FIGURE 5
The average (squared, absolute-value) Fourier spectrum (red) and the resulting ADP inter-pulse frequencies (orange) as derived from the spectrum integral (blue) for each cardiac-tissue model (rows).

Frontiers in Network Physiology
frontiersin.org

Designing a genetic algorithm
The restriction to 5-pulse pacing protocols left us with 9 degrees of freedom (5 pulse amplitudes and 4 intervals) for optimisation. We chose to use a genetic algorithm (GA) to circumvent the potential workload a brute-force approach would have entailed, which meant we had to agree on a suitable implementation for our optimisation purposes.

General structure
GAs attempt to emulate the evolutionary process in nature to improve upon known solutions to a problem with a measurable performance metric. While there is no definitive canonical GA, the most general approach works as follows (Eiben and Smith, 2015): 1. Take members from the current "population" of known solutions 2. Apply operators like mutation and/or crossover to generate new potential solutions 3. Assimilate the new solutions into the population In our case, "solution" refers to a concrete pacing protocol. Each consecutive iteration of this scheme is called a generation.
What exactly each step entails is up to a given implementation and usually tailored to a specific problem. Despite this freedom, our particular implementation had to ensure an appropriate degree of population diversity and ergodicity (i.e., the algorithm's ability to access a significant fraction of the optimisation landscape) throughout the generations: An overemphasis on population diversity would have yielded what is basically a directionless random search, and neglecting it entirely would have promoted strict elitism by favouring a restrictively narrow class of solutions instead. Either of these two cases would have compromised our GA's efficacy in finding better pacing protocols.
The price for using a heuristic lies in the abundance of parameters introduced; these are called hyperparameters to distinguish them from those pertaining to the underlying optimisation problem (i.e., the pacing protocols in the context of this study) itself. GA performance is closely tied to proper hyperparameter fine-tuning, which can be an arduous process depending on the hyperparameter count of a given GA. What follows in the next few subsections are the details of our ("the") GA, including its hyperparameters and our attempt to account for their potential influence on GA performance.

Mutation and crossover
In a given generation of the GA, there is a fixed chance (defined by the probability of mutation hyperparameter) that a mutation FIGURE 6 ADP protocol dose-response curves for our cardiac-tissue models (bottom) with a separate detailed view of the high-TR region (top). Dots mark each model's PC thresholds at TR values of 90%, 95%, 99%, and 100%, respectively. Each curve was sampled at 101 equidistant amplitude values with 10 5 termination attempts across 10 initial conditions.

Frontiers in Network Physiology
frontiersin.org operator will be consulted to generate a new protocol. When this happens, the protocol to be mutated is drawn randomly from the population and receives exactly one change in one of its parameters (Figure 7). The change in question is a replacement operation where the parameter's value is replaced by one drawn from a pre-defined interval in a uniform fashion. Each type of protocol parameter (e.g., amplitudes or recovery intervals) has its own mutation interval, and they all constitute hyperparameters as well.
This particular mutation operator may introduce drastic (in a relative sense) change to a given parameter as it employs replacement instead of scaling. However, implementing it in this way stood to ensure a consistent driving force for ergodicity and should-in theory-have allowed for the traversal of local optima at any point during the population's "evolution." This mutation operator retains full access to the entire optimisation space no matter the state of the population, which may not always be the case for scaling approaches.
Should random chance decide against mutation as a means to generating a new protocol, the GA defaults to a crossover operator instead. This operator generates two new protocols instead of one: The crossing of two random protocols begins with the selection of a random pulse number, excluding the final pulse to guarantee the results will not be two identical copies of the input protocols. Then, the two "parents" are split and recombined into two new protocols at the cutoff point defined by said pulse number (Figure 8).
"Parent" protocols are always drawn in a random, uniform fashion from the current population without replacement. This uniform approach should imply that no selection pressure is present in the choice of source protocols itself. Less selectivity in this context usually means greater ergodicity at the potential cost of convergence speed, but we could justify this choice with the highdimensional nature of the protocol-optimisation landscape (and thus inherent difficulty in its traversal) at hand.

Selection mechanism
The GA uses a percentile-based selection process to decide whether a given protocol is to be introduced to the population or discarded: By comparing the performance metrics of the candidate (i.e., its PC and TR values) to those of the existing population, it first establishes the would-be percentile rankings for both. If and only if both rankings are sufficiently high, the protocol then passes selection. The two acceptance percentile thresholds to be cleared are both hyperparameters.
An accepted protocol replaces an existing one in the population, keeping the population size (another hyperparameter) constant. The TABLE 4 ADP reference protocols. Each model has one sequence of pulse periods in-between consecutive pulses as derived from its Fourier spectrum. For each model, the protocols' uniform amplitudes were determined according to four pre-defined TR thresholds by consulting their respective dose-response curves ( Figure 6).

Model
Cutoff

FIGURE 7
Sketch of the mutation operator. A random protocol parameter (top) is selected and replaced with a randomly chosen value from a pre-defined interval (bottom). As we only considered protocols of constant pulse length, the GA implementations only ever mutated either amplitudes or recovery intervals.
Frontiers in Network Physiology frontiersin.org GA determines the candidate for replacement in a two-stage process: First, it sorts the population by one metric (PC or TR) and then by the other. This establishes a ranking where the first metric of choice serves as the principal determinant of overall protocol fitness, while the second one is effectively a tie-breaker. The protocol at the bottom of this ranking is consequently replaced by the accepted protocol. This assimilation scheme is biased by design, as it implicitly weighs the metric chosen for the first sorting more than the other. This ensures that even protocols excelling in one particular metric (potentially at the cost of the other) may still be subject to timely replacement. This is desirable behaviour as the selection process is percentile-based and accounts for both metrics, meaning outliers of this kind could reasonably skew the optimisation process. Avoiding this scenario should therefore counteract premature convergence. By alternating the order of metrics every time it consults such a ranking for replacement, the GA can apply this bias to either metric more or less equally.

Results
With a heuristic and reference protocols in hand, we ran computer simulations to study the GA's viability for uncovering protocol-performance optimisation potential across our four cardiac-tissue models of choice.

Improvements over ADP protocols
Across all four models, our top GA-optimised protocols outperformed their ADP equivalents in terms of PC efficiency only at TR values of 99% and above (Figure 9) with some GA hyperparameter sets. At 100% TR, peak GA-mediated performancegains across all hyperparameter sets were highest (Table 5), with, e.g., the FK3 ADP-protocol permitting a~80% reduction in PC through its top GA counterpart. Comparatively, the BOCF model only saw improvements of~24% over its ADP reference at peak efficacy.
Likely due to excessive selection pressure maximising the termination rate, most GA populations approached a PC-TR saturation curve that falls beneath the respective ADP doseresponse, implying overall worse PC efficiency at comparable TR. The GA's "catching up" to overall ADP-protocol performance at higher TR values may have been owed to the asymptotic nature of the ADP dose-response curves ( Figure 6): There, a TR approaching 100% demands disproportionally higher PC, enlarging the optimisation space for discovering potentially more efficient alternatives by means of GA optimisation.
In contrast, the GA achieved consistent performance gains (up to~80% PC reductions at peak TR) within the initial populations for each model. This indicates that the overall heuristic does, in fact, operate as expected in optimising a given protocol population toward more efficient pacing schemes.

Population-convergence behaviour and robustness
We identified three recurring convergence patterns in the populations' mean performance metrics, PC and TR, over the hundreds of GA evaluations conducted ( Figure 10): asymptotic, rebounding asymptotic, and quasi-linear convergence. While asymptotic convergence was the most common, its rebounding variant could be seen as evidence that the GA is not a fully "greedy" algorithm and should be capable of navigating local optima to some extent.
Metric convergence behaviour was stable under different initial populations ( Figure 11) for all models. A direct link between hyperparameters and metric convergence is difficult to establish, however, since our two population-initialisation schemes likely influenced the outcome to some degree. As we only employed one of these initialisation methods for a given model, we lack the means to compare the two and adequately quantify this potential bias after the fact.
The substantial PC and TR standard-deviations for some populations in Figures 9, 11 show that not every hyperparameter combination managed to fully converge within the allotted 10 4 generations. Fully converged data points cluster chiefly along an Sketch of the crossover operator. Given two randomly chosen source protocols (top, columns), a random pulse (the third one in this case) beyond the first is chosen as the cutoff point for the mergers of the resulting partial protocols (colour-coded). This results in two new protocols for selection (bottom).

Frontiers in Network Physiology
frontiersin.org implied saturation curve along the TR axis where the lower end of PC values lie; while this could be an indicator that the GA is principally operating in the correct direction of the PC-TR optimisation plane, the question remains whether these partially converged GA setups would have yielded superior results than those of their more rapidly converging counterparts eventually.

Emerging protocol patterns
We could not identify any overarching, reproducible protocol patterns among the GA-optimised top performers. Even reruns using the same initial population and GA hyperparameters would yield different protocol structures due to statistical variance alone ( Figure 12).
Our GA-mediated protocols likely fall into one of potentially many fine-grained local minima, past any generalisable structures of note. This could also tie in with the GA's overall incentive to optimise toward a (mathematically strictly optimal) TR of 100%, which may overemphasise overly specific protocols to the possible detriment of generalisable, similarly high-performing (e.g.,~90% TR) patterns.

Summary and discussion
With a grasp of the GA's overall performance and optimisation efficacy, we can now put our results in the larger context of our goal of protocol-performance optimisation and identify potential avenues for improvement and future research into this topic.

FIGURE 9
Overview of the GA-optimised populations for a single initial population (grey marker) per model (rows). These populations (comprising 50 protocols each) only differed in their choices of GA hyperparameters and are compared by their mean PC and TR values. Error bars indicate each metric's standard deviation. Lastly, dashed red lines show the performance of ADP protocols at a given PC value.
Frontiers in Network Physiology frontiersin.org We successfully established pacing-protocol optimisation potential between~4% and~80% over our similarly effective (in terms of TR) ADP reference protocols (Table 5). While we have managed to improve on our GA's initial populations on the order of~80% overall PC reduction at peak TR, the overall lacklustre comparison to otherwise high-performing ADP protocols (Figure 9) shows room for improvement in the general optimisation heuristic. There may also be generalisable protocol patterns to be discovered within the realm of close-to-peak-efficacy (~90% and above) protocols as a whole.
The main insight from our findings thus concerns the idea that optimisation toward a 100% success rate may be overly restrictive as far as protocol-pattern discoverability is concerned, and an equivalence (in terms of protocol fitness) among TRs abovẽ 90% may present a worthwhile avenue of investigation in the future.
The GA implementation was based on mutation and crossover operators with percentile-based selection for the assimilation of new protocols into a given population. It proved consistently capable of significantly improving all models' initial protocol-populations, but the lack of recurrently emerging protocol patterns leaves us without concrete insight into the general properties and patterns of an efficient pacing protocol.
We used only four cardiac-tissue models: Three FK and one BOCF, along with arbitrarily selected parameter sets. Each model was considered in isolation, while cross-optimisation could perhaps have helped guide the optimisation toward more general protocol structures. While we managed to replicate some of the GA's behaviour across these individual models, this was done using simplified two-dimensional simulations of fibrillation on square sheets of tissue only. Naturally, more realistic geometries along with physiological cardiac-tissue models might yield results more immediately applicable to traditional experimental settings.
Our pacing model was based on homogeneous cardiac tissue with randomly distributed, square-shaped injection sites serving as an abstraction of local tissue heterogeneities responsible for emitting wave fronts when exposed to potent electric fields as utilised by LEAP. The number of sites, their sizes, and their shape(s) were all set to rather arbitrary values for simplicity's sake, leaving yet another potential avenue for further investigation as to their effect(s) on GA performance. Perhaps more sophisticated site-distribution algorithms beyond random placement could help mimic real heterogeneity distributions more closely in this context, but the indirect implementation of virtual electrodes through the simulation 3.14 −4 100 +0

FIGURE 10
Samples for each of the three recurring metric-performance patterns we observed over the course of the GA evaluations.
Only PC examples are shown because TR convergence generally exhibited similar patterns (albeit mirrored horizontally).
Frontiers in Network Physiology frontiersin.org

FIGURE 11
Comparison of population convergence under different initial populations (columns) as well as models (rows). Each final population consisted of 50 pacing protocols and is placed according to its mean PC and TR values. Error bars indicate the corresponding standard deviations.

FIGURE 12
Five examples of top GA-optimised protocols for two FK models (rows). These protocols had the highest TR with the lowest PC among their respective populations. Said populations were generated by a single 10 4 -generation application of our GA using one set of hyperparameters and initial population per model.
Frontiers in Network Physiology frontiersin.org 14 of actual heterogeneous tissue domains may also be worth considering.
Lastly, the particular GA setup we employed is merely one of many possible alternatives, and others may prove yet more effective at optimising pacing protocols. We used either mutation or crossover to generate one or two new protocols in a given generation and selected the respective source protocol(s) entirely randomly every time, which deviates from the more "canonical" approach (phrased in the context of pacing protocols) of using, e.g., tournament selection (Blickle and Thiele, 1996) to evolve the protocol population.
The next step in computer-aided pacing optimisation could entail the cross-optimisation with multiple models, fitness equivalence of TR values beyond~90%, transition to threedimensional simulations, more sophisticated pacing models, the testing of alternative GA implementations, or a combination thereof. We established a foothold in pacing optimisation through the outline of this paper, but more work into studying the GA, its behaviour, and its shortcomings is needed before considerable time and effort is spent working with, e.g., physiological models in more complex geometries beyond square sheets of cardiac tissue.
Other popular heuristics like simulated annealing (Ingber, 1993) may also warrant investigation for their efficacy in protocol optimisation based on our pacing model. It would at the very least be pertinent to compare their results to our GA's, especially in light of the convergence issues we have uncovered over the course of the GA-based optimisation effort.
Further optimisation of anti-fibrillation control may be achieved by proper timing of the application of the pulse sequence, as it has recently been shown that the timing of single-shock termination has a significant influence on the probability of success (Ji et al., 2017;Steyer et al., 2023).
We can conclude this paper with a successful foray into pacingprotocol optimisation where we showcased the optimisation capabilities of a simple genetic algorithm for randomised protocol populations, as well as substantial pacing-energy reductions over ADP protocols at 100% sucess rate. In the process, we identified multiple avenues of investigation for future projects; said projects could focus on general improvements to the GA design for better protocols or taking the simulations to a more realistic level to better facilitate their translation into experiments on live cardiac tissue for verification purposes.

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

Author contributions
MA, UP, TL, and SL conceived the study, MA performed the simulations and analyzed the data. MA, UP, TL, and SL wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding
The authors acknowledge funding by the Excellence Cluster MBExC (project "Motion analysis of excitable biological media using machine learning"), by the Max Planck Society (MPG) and the German Centre for Cardiovascular Research (DZHK) e.V.