A Reproducible Protocol to Assess Arrhythmia Vulnerability in silico: Pacing at the End of the Effective Refractory Period

In both clinical and computational studies, different pacing protocols are used to induce arrhythmia and non-inducibility is often considered as the endpoint of treatment. The need for a standardized methodology is urgent since the choice of the protocol used to induce arrhythmia could lead to contrasting results, e.g., in assessing atrial fibrillation (AF) vulnerabilty. Therefore, we propose a novel method—pacing at the end of the effective refractory period (PEERP)—and compare it to state-of-the-art protocols, such as phase singularity distribution (PSD) and rapid pacing (RP) in a computational study. All methods were tested by pacing from evenly distributed endocardial points at 1 cm inter-point distance in two bi-atrial geometries. Seven different atrial models were implemented: five cases without specific AF-induced remodeling but with decreasing global conduction velocity and two persistent AF cases with an increasing amount of fibrosis resembling different substrate remodeling stages. Compared with PSD and RP, PEERP induced a larger variety of arrhythmia complexity requiring, on average, only 2.7 extra-stimuli and 3 s of simulation time to initiate reentry. Moreover, PEERP and PSD were the protocols which unveiled a larger number of areas vulnerable to sustain stable long living reentries compared to RP. Finally, PEERP can foster standardization and reproducibility, since, in contrast to the other protocols, it is a parameter-free method. Furthermore, we discuss its clinical applicability. We conclude that the choice of the inducing protocol has an influence on both initiation and maintenance of AF and we propose and provide PEERP as a reproducible method to assess arrhythmia vulnerability.


INTRODUCTION
Atrial fibrillation (AF) is the most frequent cardiac arrhythmia and a progressive pathology associated with high morbidity and mortality (Hindricks et al., 2020). Despite recent advances in both diagnostic and therapeutic techniques, the success rate for the standard-of-care treatment, catheter ablation, is sub-optimal in patients with persistent AF (Verma et al., 2015;Morady and Latchamsetty, 2018). The modest efficacy reflects the complexity of the underlying phenomena and our incomplete understanding of the mechanisms of initiation, maintenance and progression of AF episodes (Andrade et al., 2014;Mann et al., 2018). In clinical practice, electrical stimulation has been widely used to diagnose and guide therapy of arrhythmias. Despite the high prevalence of atrial rhythm disorders, the sensitivity, specificity and reproducibility of mostly ventricular stimulation has been investigated (DiCarlo et al., 1985;Kudenchuk et al., 1986). An equally critical evaluation of protocols for induction of AF and atrial flutter is currently lacking. Kumar et al. (2012) highlighted the importance of the choice of the protocol used to induce AF even in patients without AF history and/or structural heart disease. The incidence of initiation and maintenance of AF varied according to gender, method of induction and number of inductions. In addition, the stimulation sites, pacing methods, number of AF inductions, use of pharmacological provocation and the definition of inducibility based on AF duration vary among different studies (Hakan et al., 2004;Essebag et al., 2005;Jaïs et al., 2006;Richter et al., 2006). Currently, various protocols are used to test AF inducibility before and after ablation procedures. For instance, some groups have used burst pacing at a fixed cycle length, while others have paced at the shortest cycle length which resulted in loss of 1:1 capture. Kumar et al. (2012) observed that the incidence of inducible or sustained AF was significantly higher with decremental pacing compared to burst pacing, as was the total duration of induced AF. They concluded that the adoption of AF inducibility as final electrophysiological endpoint is critically dependent on the variations in the definition of inducibility, aggressiveness of AF induction protocol and the number of AF inductions. However, they could not compare burst vs. decremental pacing within the same patient as this is not feasible clinically for ethical reasons. In another clinical study (Krol et al., 1999), a programmed atrial stimulation protocol for induction of sustained arrhythmia was evaluated. They were able to induce and maintain arrhythmic episodes using a train of pulses close to the effective refractory period in 39/44 (89%) patients. The study demonstrated that employing only two atrial sites and three atrial extra-stimuli induced either AF or atrial flutter in 89% of the patients with previous AF history and in 7% of the control group without documented arrhythmias. Moreover, inducibility was proposed as a predictor of long-term AF recurrence (Essebag et al., 2005).
Computational modeling has been proven to be a useful tool for assessing arrhythmia vulnerability (Arevalo et al., 2016;Zahid et al., 2016;Azzolin et al., 2020) and for supporting ablation planning (Lim et al., 2017;Boyle et al., 2019;Loewe et al., 2019;Roney et al., 2020). However, different protocols used to induce arrhythmia in simulations (Krummen et al., 2012;Matene and Jacquemet, 2012;Bayer et al., 2016;Roney et al., 2016Roney et al., , 2018Zahid et al., 2016;Boyle et al., 2019) are not only making studies difficult to compare, but are also influencing the decision on whether an atrial model is vulnerable to AF and are therefore crucial for identifying the optimal ablation targets. In Boyle et al. (2019), arrhythmic episodes were induced by pacing from 40 evenly distributed sites applying a train of 12 electrical stimuli with decreasing cycle length from 300 to 150 ms. Their choice of protocol was motivated by a previous publication , in which, however, AF was initiated by delivering five triggering ectopic beats from the right superior pulmonary vein with a fixed coupling interval of 400 ms with the sinus rhythm but variable basic cycle length of 155 or 160 ms, depending on inducibility. In Bayer et al. (2016), each pulmonary vein was individually paced with five beats at a fixed cycle length but different coupling intervals. The values for coupling interval and cycle length were according to Haissaguerre et al. (1998). Roney et al. (2018) tested arrhythmia inducibility via extra stimuli coming from each pulmonary vein. They already highlighted the influence of variable cycle length, pacing intervals and pacing location on arrhythmia initiation. Conversely, Matene and Jacquemet (2012) proposed to induce AF by manually placing 1-6 phase singularities on the atrial surface, reconstructing an activation time map using an interpolation algorithm based on the eikonal equation and using that as initial condition for a monodomain simulation.
Considering the increasing use of catheter ablation and pacing therapies for atrial arrhythmias, which often take the suppression of inducibility as one of the endpoints of treatment, the definition of a commonly acknowledged electrical stimulation protocol is needed. A standardized methodology is critical since the use of different protocols could lead to contrasting conclusions on whether or not a patient or a patient-specific digital twin model is vulnerable to arrhythmia. In this study, we quantitatively evaluate different methods to induce arrhythmia, investigate their impact on both initiation and maintenance of AF episodes and propose an easily reproducible protocol to assess AF vulnerability of a particular atrial model.

Atrial Models
In this work, two highly detailed volumetric bi-atrial geometry derived from magnetic resonance images were used (Krueger et al., 2013). The tetrahedral mesh had an average edge length of 0.5 mm. Fiber orientation was calculated by a semiautomatic rule-based algorithm (Wachter et al., 2015). We updated the models used in Krueger et al. (2013) by including a pectinate muscle running from the sinus node to the apex of the right atrial appendage to provide a fast activation of the appendage. Different conductivity and anisotropy values were implemented in working myocardium, pectinate muscles, bachmann bundle, inferior isthmus and crista terminalis to consider the heterogeneity in the atria (Loewe et al., , 2016. Conductivities were tuned to reach four different total activation times in sinus rhythm of 130, 151, 179, and 199 ms and respective global average conduction velocities of 0.7, 0.66, 0.55, and 0.49 m/s, respectively. We will refer to these models as H1, H2, H3, and H4. In the second geometry, we applied the same mean global conduction velocity as in model H4. The resulting atrial model had a total activation time of 152 ms and will be referred to as H4B. The total activation times were chosen to be all ≥130 ms since patients with a P-wave duration longer than 130 ms were shown to have higher risk for AF (Lemery et al., 2007;Nielsen et al., 2015). Depending on the simulated scenario, the original Courtemanche et al. (1998) model (H1, H2, H3, H4, and H4B) or a variant reflecting AF-induced remodeling (UII and UIV) (Loewe et al., 2014) represented the myocyte membrane dynamics. Fibrotic tissue was included in the regions of both atria most-often exhibiting fibrotic substrate in patients with AF, which are pulmonary veins antrum, left lateral wall, anterior wall, posterior wall and septum in the left atrium (Akoum et al., 2012;Benito et al., 2018;Higuchi et al., 2018) as well as right atrial appendage and septum in the right atrium (Cao et al., 2010;Akoum et al., 2012).
Two different fibrosis severity stages were implemented: Figure 1A shows the Utah stage II case (UII), in which 19% of the left atrial wall (LAW) and 5% of the right atrial wall (RAW) were modeled as fibrotic. In the most severe case (Figure 1B), 39% of the LAW and 11% of the RAW was modeled as fibrotic, classifying as Utah stage IV (UIV). To account for structural remodeling and the presence of scar tissue, we set 50% of the elements in the fibrotic regions as almost not conductive (conductivity of 10 −7 S/m). This approach modeled the macroscopic passive barrier behavior caused by the electrical decoupling of the myocytes in the tissue infiltrated by fibrosis, also referred to as "percolation" . In the other 50%, several ionic conductances were rescaled to consider effects of cytokine-related remodeling  (−50% g K1 , −40% g Na , and −50% g CaL ).

Protocols to Test Inducibility
We tested inducibility of arrhythmic episodes by pacing from stimulus locations placed on the endocardial surface with interpoint distance of 1 cm, resulting in up to 227 evenly distributed points. We systematically evaluated and compared three stateof-the-art methods: phase singularity distribution (PSD), rapid pacing (RP) and pacing at the end of the effective refractory period (PEERP). In all pacing protocols (RP and PEERP), a transmembrane current of 30 µA/cm 2 was injected in an area of 2 × 2 mm. We considered a point to be "inducing" if the application of one protocol at that location, induced and sustained an arrhythmia for at least 1.5 s after the end of the protocol. All the protocols used in this study are available open source in the examples section of openCARP (Vigmond et al., 2003;Sánchez et al., 2020;Augustin et al., 2021) at www.opencarp.org/documentation/examples.

Phase Singularity Distribution
The phase singularity distribution (PSD) method places phase singularities (PSs) in the atria, constructs an activation time map by solving the Eikonal equation and finally uses this as initial state for a monodomain excitation propagation simulation (Jacquemet, 2010;Matene and Jacquemet, 2012). The parameters of this method are the cycle length of the reentrant wave, the number of PSs to initialize and the rotation direction of each PS. We set a single PS rotating in anti-clockwise direction (looking from the endo-to the epicardium) in one of the 227 points for each of the 227 simulated scenarios. We seeded only one PS for the sake of consistency with the pacing protocols in which we stimulated from one location per simulation. The cycle length was set to 315 and 168 ms, for the original (Courtemanche et al., 1998) model and the chronic AF variant (Loewe et al., 2014), respectively (5% longer than the effective refractory period). However, this method can be applied only in simulations.

Rapid Pacing
The rapid pacing (RP) protocol consists of a train of stimulation pulses with decreasing coupling intervals (CI). We extended the protocol by supporting a variable number of pulses with the same CI before the next decrement. Moreover, we implemented the option to check for successful initiation of an arrhythmia after every stimulation instead of only at the end of the protocol. This extended RP protocol is defined by RP N, {B,E} s,l where N is the number of stimuli with the same CI, {B, E} defines if we were checking for the induction of an arrhythmia after every beat (B) or only at the end of the protocol (E), s is the starting CI and l is the last CI, both in ms. The CI was continuously decremented from s to l in steps of 10 ms. s was different for each action potential phenotype: 300 ms for the control case and 200 ms in the AF-induced remodeling setup. l was also distinct between control and AF remodeled and chosen as 200 and 130 ms, respectively. N was incremented from 1 to 4. For example, the RP protocol with N = 1, s = 200 ms, l = 130 ms and arrhythmia checking only at the end resulted in the following train of pulses: which is the most common RP protocol (Krummen et al., 2012;Zahid et al., 2016;Boyle et al., 2019). The case N = 2 was as follows: 200, 190, 190, . . . , 140, 140, 130, 130 check}.
(2) The RP protocol with N = 1, s = 200 ms, l = 130 ms and arrhythmia checking after every beat yields: (3) To summarize, the RP protocols were tested with N increasing from 1 to 4 and coupling interval between 300 and 200 ms for the control case and 200 and 130 ms for the AF remodeled case. The parameters of the RP protocols are the maximum number of beats with the same CI, first and last CI, the CI decrement, and the stopping criteria (arrhythmia checking in between or at the end of the protocol). Moreover, the RP is currently used in clinical practice.

Pacing at the End of the Effective Refractory Period
The pacing at the end of the effective refractory period (PEERP) protocol triggers stimuli at the end of the effective refractory period. Therefore, each stimulus was delivered as soon as the underlying tissue had recovered from the previous activation and was able to initiate a new wave propagation. We implemented a run-time binary search method to find the minimum time at which a new depolarization wave could locally spread (effective refractory period) with a temporal resolution of 1 ms. We tested the influence of a temporal ERP uncertainty up to 8 ms in the model UIV. Based on initial experiments, the atrial action potential duration at 94% repolarization computed during sinus rhythm was chosen as first guess of the effective refractory Black regions indicate fibrotic tissue in which 50% of the elements were non-conductive and the remaining 50% electrophysiologically remodeled due to cytokine effects.
period. The criterion for successful wave propagation was a transmembrane potential ≥ −50 mV in at least one node in a ring 4-6 mm around the stimulation site. In this study, the maximum number of beats delivered from one stimulation site was set to 4. Consequently, the stimulation at one pacing location was stopped if the maximum number of beats was reached or if an arrhythmia was initiated. The only parameter for this new PEERP protocol is therefore the maximum number of stimuli to apply at a specific location. We will discuss in section 4.4 how the PEERP could be implemented during in-vivo procedures.

Arrhythmia Classification
Reentrant dynamics can, in many cases, be simplified by considering the behavior of their organizing centers-the PSs. We followed the method presented by Clayton et al. (2006) to detect PSs and track the spatio-temporal behavior of each PS during the whole simulated episode. We defined a PS as longliving PS (llPS) if it was sustained for at least 500 ms. Moreover, we distinguished the llPS which were stable within an area enclosed in a bounding box of 5 cm edge length (stable llPS) for the last simulated 1.5 s, from the ones which were meandering largely (non-stable llPS). The episodes in which no llPS was found were classified in multiple wave fronts (Multi) or flutter (Fl) by visual inspection. We, therefore, classified each arrhythmic episode into four categories: Multi, Fl, non-stable llPS, and stable llPS. The Multi class was defined as a complex and disorganized series of multiple wave fronts merging and colliding with each other without a clearly identifiable llPS. Periodic macro-reentries were interpreted as Fl cases. An example of stable llPSs and multiple wave fronts episodes are shown in Figure 2.

Subdivision of Atria in Spatial Segments
To localize the inducing points and drivers sustaining atrial activity, we partitioned the atria into different regions. The atria were subdivided into 28 segments (Figure 3): 19 in the left atrium [four pulmonary vein segments (LIPV, RIPV, LSPV, and RSPV), four segments at the roof, five segments on the posterior wall, two segments in the septum, four segments on the anterior wall] and nine in the right atrium (inferior and superior vena cava rings, coronary sinus, cavotricuspid isthmus and septum, sinus node, right atrial appendage, anterior wall, four segments on the posterior wall of the RA).

Computational Tools
The spread of the electrical activation in the atrial myocardium was simulated by solving the monodomain system using openCARP (Vigmond et al., 2003;Sánchez et al., 2020) and a time step of 0.02 ms.

Inducibility
The number of inducing points and the relative frequency of arrhythmic mechanisms induced in each model are shown in and RP E 200,130 for the models UII and UIV. However, in Figure 4 they are summarized as RP B and RP E . The number of inducing points increased with both a lower average global conduction velocity and higher amount of fibrotic tissue. No arrhythmic episode could be induced and sustained in model H1, which is why this model is not included in the results. PEERP was the only method which induced all the possible arrhythmic mechanisms described in section 2.3 in most of the models (4/5). The application of the RP protocol resulted in inducing a majority of atrial flutter mechanisms in the model UII and no multi-frontal episode in either model with fibrosis. Moreover, the PEERP method initiated arrhythmia from the most segments, averaged over all models, as displayed in Figure 5. However, RP B was the only method which could induce and sustain arrhythmia pacing from all the segments in the UIV model. Nonetheless, inducing points belonging to only 4/28 segments were found applying the RP protocol in the model UII. The atrial segments in which inducing points were found to initiate and maintain stable llPSs applying the different protocol can be found in the Supplementary Material. In the left atrium, the majority of inducing points were found in segments 4,5,7,9,18,19. In the right atrium, the segments 25, 26, 27 were the most inducible.
The CI of RP B at which arrhythmia was initiated was around 130 and 250 ms for the models with and without structural remodeling, respectively ( Figure 6A). The following tendency was observed: a slight increase of the inducing CI with lower global conduction velocity and higher amount of fibrosis. The inducing number of stimuli with the same CI (N) was in most of the cases N = 2 for both RP B and RP E , as presented in Figure 6B. However, in the model UIV N = 1 was often sufficient to initiate arrhythmic episodes using RP B . The required total number of stimuli to induce arrhythmia by each protocol is displayed in Figure 6C. The methods PEERP, RP B and RP E required on average 2.7, 13.6, and 17.6 beats respectively to initiate an arrhythmic episode. The different protocols needed a similar total number of beats when applied to the model H4B, as shown in the Supplementary Material.
Furthermore, the application of PEERP resulted in an increasing number of segments in which inducing points  were found with both decreased global CV and greater amount of fibrotic tissue, as shown in Figure 5. On one hand, the RP protocols were the only ones to initiate and sustain arrhythmia pacing from almost all the segments in the model UIV. On the other hand, both RP B and RP E could induce pacing from only four segments in the model UII.
Moreover, the simulation time required to complete each protocol for one stimulation point was 3.0, 19.5, 8.4 s for PEERP, RP B , and RP E , respectively (Figure 7).
The number of inducing points found in the model UIV remained stable between 141 and 146 out of 227 stimulation points when increasing the time tolerance for the local effective refractory period computation up to 8 ms as shown in Figure 8.
Additionally, subsets of the initial pacing locations were extracted to evaluate the sensitivity of each protocol to a reduction of the number of pacing locations, i.e., an increase of the inter-point distance. The sensitivity was computed as the percentage of stable llPSs found when pacing from the initial pacing points with 1 cm inter-point distance. We  considered increasing inter-point distances of 1.5, 2, 2.5, and 3 cm. The sensitivity analysis results are presented in Table 1. When applying the protocols only at pacing points located 3 cm away from each other, leading to a total of only 11 bi-atrial pacing sites, all the protocols could reproduce more than the 90% of the locations maintaining stable llPSs in the model UIV. The PEERP and the RP B were the only methods to identify some of the stable llPSs sustaining areas with a greater inter-point distance when applied to all models.

Maintenance
All protocols showed a tendency for increasing number of segments in which stable llPSs were sustained in models with lower global CV or higher percentage of fibrotic remodeling (Figure 9). PEERP was the protocol which sustained most stable llPSs in the models without structural remodeling. In the UII and UIV models, the PSD method was the one identifying more segments vulnerable to maintain stable llPSs, followed by PEERP. The RP protocols classified subsets of atrial segments found by the other methods as vulnerable to sustain stable llPSs in most of the cases.
The PEERP method was the only one identifying the left atrial roof (segment 9) as vulnerable for maintenance of stable llPSs in all the models. All the pacing protocols (PEERP and RP) showed most of the stable llPSs in segments located in the left atrium. The inclusion of fibrotic remodeling (models UII and UIV) led to the maintenance of stable llPSs in the segments containing fibrosis. However, PSD classified as vulnerable to maintain many segments located in the right atrium too. Moreover, the PSD method yielded many stable llPSs sustaining in segments without fibrosis in the models UII and UIV. Independent of the choice of the inducibility protocol, most of the stable llPSs were maintained in the left posterior wall, roof, and septum (segments 5, 6, 9, 14, 16, 25, 27). The atrial segments in which stable llPSs were sustained after applying the different protocols can be found in the Supplementary Material. The distance between the centers of the path traveled by stable llPSs and their corresponding inducing points was around 5-7 cm in most models using PEERP and RP B (Figure 10). Applying PSD and RP E the distance was between 2.5 and 5 cm in the models without AF specific remodeling and between 5 and 7 cm in the others.
All induction protocols sustained on average one stable llPSs per simulation, as shown in Figure 11.

Inducibility
Each of the four protocols was applied at 227 locations for each of the seven atrial models, resulting in 6,356 3D bi-atrial electrophysiological simulations. This is, to the best of our knowledge, one of the most extensive computational studies on the influence of pacing protocols on AF inducibility and maintenance. Regarding the vulnerability to induce arrhythmia, we demonstrated that the PEERP method provokes a bigger variety of different mechanisms compared to the other protocols. PEERP was able to induce diverse degrees of arrhythmia complexity, ranging from a flutter mechanism to chaotic multiple wave fronts. PEERP unveiled these complex mechanisms in most of the models which could potentially sustain them, even in the cases in which it was not the protocol with the highest number of inducing points. Furthermore, we showed that PEERP yields, overall, the highest number of segments from which arrhythmia could be induced (Figure 5), meaning that this protocol is less dependent on the pacing location.
The results of the sensitivity analysis (Table 1), confirmed the ability of the methods PEERP to identify a good percentage of the total stable llPSs in all the models even with pacing locations placed at a greater inter-point distance between each other. These results are in line with Krol et al. (1999), who could induce arrhythmia in 89% of patients with AF history with a train of stimuli close to the effective refractory period from a couple of pacing locations.
Moreover, we tested the influence of a temporal uncertainty up to 8 ms for the effective refractory period on the final number of inducing points (Figure 8). Since the number of inducing points remained relatively constant when increasing the temporal tolerance, the PEERP method can be used to assess atrial model vulnerability to AF with a temporal resolution up to 8 ms without affecting the results markedly.
The application of the PEERP method showed that a similar number of beats (2-3) was needed to initiate arrhythmia in all the models used in this study. This result matches with the findings by Krol et al. (1999), where three extra-stimuli were enough to initiate either AF or atrial flutter in patients both with and without history of AF. When fixing the maximum number of beats to 3, PEERP becomes a parameter-free method suitable to test inducibility in various electrophysiological models. Moreover, the smaller simulation time needed to perform the entire PEERP method from one stimulation point could give the possibility to test inducibility from more pacing locations to ensure the identification of all possible arrhythmic episodes which can arise in a given atrial model.
Moreover, we showed a strong dependence of the induced episodes on the CI of the RP protocol (Figure 6). Moreover, this FIGURE 9 | Number of maintaining atrial segments after applying each protocol in all the 227 stimulus locations. We simulated 1.5 s after the end of each protocol to observe the evolution of the arrhythmia. parameter had to be tuned for each ionic model. We noticed that a smaller CI was needed to induce arrhythmia in models less vulnerable to initiate AF (higher CV or limited fibrosis). This is in agreement with the need of being more aggressive to induce AF in patients belonging to less severe AF subgroups. The small number of llPSs initiated by the RP in the model UII clearly showed the limitation of user defined CI window. Furthermore, we noticed that a single beat per CI is usually not enough to induce arrhythmia and then a second beat with the same CI is needed. More than two beats with the same CI rarely induced arrhythmia. This could confirm the clinical finding of Kumar et al. (2012), who found a protocol with decreasing CI (decremental pacing) to be more effective in inducing AF than a fast train of pulses with fixed CI (burst pacing).
Our study highlighted the importance of checking for arrhythmia initiation after every beat. Without this as a stopping criterion for the protocol, many arrhythmic episodes induced by previous stimuli were terminated by subsequent stimuli. Indeed, RP B initiated arrhythmia before reaching the lowest CI and with a lower number of beats per CI than the predefined maximum N = 4 (Figure 6). Finally, arrhythmic episodes were always more easily initiated by applying RP B compared with RP E , however at a much higher computational cost. Our results are in line with clinical findings by Kumar et al. (2012), where inducibility differed depending on the method of induction, the number of induction attempts and the patient (in our case represented by different atrial models).
The segments in which most of the inducing points were found were also the ones including most anatomical and/or electrophysiological heterogeneity (e.g., pulmonary veins, mitral valve and posterior wall in the left atrium, crista terminalis and pectinate muscles in the right atrium).

Maintenance
The areas identified as vulnerable to sustain reentries were different between the various protocols, confirming the importance of the choice of the induction protocol not only in initiation but also in maintenance of AF. Krol et al. (1999) were able to induce AF by applying stimuli close to the effective refractory period even in patients without previous AF history. In line with these findings, PEERP-induced stable llPSs were maintained in the models without fibrotic remodeling but slower global conduction velocity.
There was a strong correlation between maintenance of stable llPSs and presence of fibrotic tissue, mostly when applying pacing protocols (PEERP and RP). As shown in Azzolin et al. (2020), the regions of the crista terminalis and pectinate muscles (located in segments 25 and 27) were prone to sustain AF in the right atrium due to their high degree of heterogeneity in fiber architecture and conduction velocity.
The PSD method labeled areas in the right atrium (segments 24 and 27) as highly vulnerable to sustain llPSs, which almost no other protocol found. This was because PSD is by definition setting PSs as initial conditions and could lead to PSs maintained close to the point in which they were initiated, mostly if not attracted by fibrotic tissue patches.
The arrhythmia check as stopping criteria played an important role not only in the initiation but also in the maintenance of stable llPSs, since it led to the identification of different segments as vulnerable to sustain reentry when comparing RP B and RP E .
The average geodesic distance of 5-7 cm between rotor cores and corresponding inducing points showed that, in most cases, stable llPSs sustain in areas not in the pacing location's neighborhood. Consequently, we needed to check arrhythmia initiation in the whole atria, not only in the proximity of the stimulation point.
The PEERP recognized the areas which we expected to be more vulnerable to sustain rotors, highly heterogeneous regions or containing fibrotic tissue. Moreover, the PEERP identified most of the segments classified as vulnerable to maintain stable llPSs by the other protocols too. This showed a high sensitivity and specificity of the PEERP method in discriminating between areas which are prone to sustain stable rotors.

Clinically Important Observations
We noticed stable llPSs dominating in atrial models including fibrotic tissue and perpetuating mostly in segments containing fibrosis. This confirmed the link between re-entrant drivers dynamics and the fibrotic tissue distribution (McDowell et al., 2015;Zahid et al., 2016). On the contrary, non-stable llPSs prevailed in atria without structural remodeling.
Moreover, we showed that not only fibrotic tissue distribution and conduction velocity, but also the induction protocol are influencing both initiation and progression of AF episodes, confirming what was clinically displayed by Kumar et al. (2012).
Different arrhythmic mechanisms were induced in our atrial models, with various degree of complexity. This could support that the rivaling theories of rotors (Krummen et al., 2015) and multiple wavelet (Moe, 1962) were both right and can co-exist. Sometimes one mechanism dominates and sometimes the other.

Clinical Applicability of the PEERP Protocol
Rapid decremental pacing and burst pacing are considered as state-of-the-art protocols to induce arrhythmia in the atria (Essebag et al., 2005;Kumar et al., 2012). Moreover, suppression of inducibility with these protocol has been used as endpoint of ablation treatment by various groups (Hakan et al., 2004;Essebag et al., 2005;Jaïs et al., 2006). However, the strong influence of the chosen protocol to initiate and maintain AF has been previously shown (Kumar et al., 2012) and was confirmed mechanistically in our work. Furthermore, the lack of a consensus regarding which method to use to test inducibility makes studies hard or impossible to reproduce and to compare. Patient hearts are electrophysiologically diverse and human atrial myocardium is intrinsically heterogeneous, calling for a pacing protocol which is as much as possible independent from human-defined parameters (e.g., basic cycle length, coupling interval, and decreasing step). A programmed atrial stimulation protocol (Krol et al., 1999) pacing close to the effective refractory period has shown positive predictive accuracy in inducing AF of 95% using only a few pacing locations and three atrial stimuli. We showed that an automatically adjustable pacing protocol stimulating at the end of the effective refractory period is able to initiate arrhythmia episodes with on average only 2-3 stimulations, in accordance with what was observed in Krol et al. (1999). However, effective refractory period is clinically normally determined by a pacing protocol, like S1-S2 or rampdown. Therefore, it is dependent on the exact pacing protocol and location. Moreover, the determination of local effective refractory period with pacing protocols will take some time and is only feasible with limited spatial resolution. In contrast, RP protocols require no setup, and in fact, may be clinically faster than local estimation of effective refractory period. However, Verrier et al. (2016) presented a method to assess atrial repolarization without provocative electrical stimuli using standard clinical catheters. This opens the possibility of applying our proposed PEERP method in clinical practice under the assumption that repolarization time can be measured accurately enough to obtain unidirectional block using, e.g., the approach suggested by Verrier et al. (2016). We showed that a temporal resolution up to 8 ms should be sufficient to assess AF inducibility.
To the best of our knowledge, the PSD method is not clinically applicable. In contrast, pacing protocols are common in clinical practice during invasive electrophysiological procedures using a pacing catheter. However, other techniques could be used to deliver electrical stimuli on demand. Crocini et al. (2016) and Scardigli et al. (2018) presented optogenetics as a feasible way to manipulate cardiac wave propagation and to discharge customized stimulation patterns across the whole heart with reaction time within 2 ms.

Limitations
Our atrial geometry did not include heterogeneous atrial wall thickness, which has been shown to influence the dynamics and stability of reentrant drivers (Azzolin et al., 2020). However, we are confident that our results hold in presence of heterogeneity in atrial wall thickness, since the shown consistency of results with the PEERP on the variety of cases we investigated also translates to new variations. In our study, only two different atrial anatomies were used. A deeper analysis on the influence of geometrical variability on AF vulnerabilty could be addressed by the use of atrial shape models (Nagel et al., 2021). We carried out monodomain simulations which could have affected the dynamics of AF, even if the bidomain equations did not show significant difference in reproducing the wave propagation in thin-walled atrial tissue (Potse et al., 2006). Due to the computational cost of this extensive study, we decided to limit the simulation time to 1.5 s after initiation of an arrhythmia. A meandering PS observed in our study could potentially stabilize later and affect the areas of maintenance. Clinically, testing of inducibility is mostly performed after ablation. The inclusion of ablation lines could potentially change the results of both inducibility and maintenance of AF episodes. They will, on one hand, limit the available space for AF to be sustained. On the other hand, different ablation patterns could provide additional substrate for reentries to be initiated and maintained. The scope of this work is to show how the choice of the induction protocol affects AF vulnerability in different atrial models. The next step could be to provide patient-specific ablation strategies to terminate the sustained AF episodes.

CONCLUSION
Our study highlights the influence of different arrhythmia induction protocols for the assessment of both AF initiation and maintenance of AF. Our newly proposed PEERP protocol offers a reproducible, comprehensive and computationally fast method to assess vulnerability. PEERP was able to provoke different degrees of arrhythmia complexity and unveil areas prone to maintain AF with a low number of stimuli, thus computationally inexpensive. The open source availability will facilitate adoption of the parameter-free PEERP method as a community standard. Therefore, we suggest the PEERP protocol (with N = 4) as a default protocol for future studies. This work is a basis to increase comparability and reproducibility of in silico arrhythmia vulnerability studies and could prove feasible to be applied clinically as well.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: www.opencarp.org/ documentation/examples.

AUTHOR CONTRIBUTIONS
LA and AL conceived and designed the study. LA constructed the atrial models, implemented the protocols, ran the simulations, analyzed the data, and drafted the manuscript. SS developed the method to subdivide the atria into segments and to compute the geodesic distance. All authors edited and approved the manuscript.