Mechanisms of arrhythmogenesis related to calcium-driven alternans in a model of human atrial fibrillation

The occurrence of atrial fibrillation (AF) is associated with progressive changes in the calcium handling system of atrial myocytes. Calcium cycling instability has been implicated as an underlying mechanism of electrical alternans observed in patients who experience AF. However, the extent to which calcium-induced alternation of electrical activity in the atria contributes to arrhythmogenesis is unknown. In this study, we investigated the effects of calcium-driven alternans (CDA) on arrhythmia susceptibility in a biophysically detailed, 3D computer model of the human atria representing electrical and structural remodeling secondary to chronic AF. We found that elevated propensity to CDA rendered the atria vulnerable to ectopy-induced arrhythmia. It also increased the complexity and persistence of arrhythmias induced by fast pacing, with unstable scroll waves meandering and frequently breaking up to produce multiple wavelets. Our results suggest that calcium-induced electrical instability may increase arrhythmia vulnerability and promote increasing disorganization of arrhythmias in the chronic AF-remodeled atria, thus playing an important role in the progression of the disease.

reentry over a narrow range of CIs, which did not sustain after reentry propagated once around the pulmonary vein (referred to as one reentrant cycle).
In contrast, all eight ectopic beat locations produced post-S2 reentrant activity in ALT slow (Fig. 3b). After the odd beat, any post-S2 activity in ALT slow terminated after one reentrant cycle, as in the ALT fast simulations. After the even beat, however, >1 reentrant cycles were induced at two S2 stimulus locations in ALT slow (RS0 and RS1). These stimuli were located adjacent to areas of early repolarization during the previous S1 beat (purple areas in Fig. 2b, bottom left). Because of the steep repolarization gradients during the even S1 beat, stimuli such as RS3, which experienced similar repolarization times but were not adjacent to areas of early repolarization, had dramatically different outcomes. With RS3, the S2 stimulus was surrounded by refractory tissue at all CIs and resulted in propagation failure.
In some cases, reentry around the pulmonary vein degenerated into fibrillation-like conduction. At RS0, three CIs resulted in reentry that broke up into multiple wavelets after about three cycles. In one simulation (RS0, CI = 270 ms), arrhythmia sustained > 2 s after the S2 stimulus, during which the initial reentry around the RSPV eventually terminated and an unstable scroll wave later appeared and meandered in the RA. These results demonstrate that repolarization heterogeneity induced by CDA resulted in increased atrial vulnerability to arrhythmia induction by ectopic beats. Arrhythmia complexity during fast pacing. We next investigated the effect of CDA on arrhythmia complexity in the ALT fast and ALT slow atria models in a new set of simulations in which wave breakup and reentry were induced during a dynamic pacing protocol. The atria were paced from the SAN region at progressively faster rates (from 750-ms to 250-ms CL, 32 cycles at each pacing rate, see Methods). Complexity was quantified by tracking filaments, which are the organizing centers of scroll waves where phase singularities (PSs) occur. Discordant APD alternans developed at 270-and 350-ms CL pacing in ALT fast and ALT slow , respectively. During 270-ms CL pacing, several discordant regions were present in both models, but APD alternans magnitude was noticeably larger in ALT slow (red and blue regions, Fig. 4a). As pacing rate increased to 260-ms CL, this repolarization heterogeneity in ALT slow led to wavebreak and filament formation near nodal lines, where regions with alternans of opposite polarity meet ( Fig. 4b and Supplementary Video S3). Filaments appeared earlier (during 270-ms CL pacing) and were present longer in the ALT slow model (Fig. 5a). During arrhythmia (i.e. at time points when filaments were present), ALT slow tended to have significantly more filaments than ALT fast (median was 2 and 6 filaments for ALT fast and ALT slow , respectively), indicating that CDA increased the complexity of arrhythmia ( Fig. 5b and Table 1).
We investigated two possible contributors to the increased number of filaments in ALT slow : additional new filament formation or increased filament duration. Figure 5c compares new filament formation during fast pacing (270-ms to 250-ms CL) in ALT fast and ALT slow . New filaments appeared steadily during 260-ms and 250-ms CL pacing in ALT slow but slowed in ALT fast after about 10 beats. Thus, to a large degree, new filament formation accounted for the increased number of filaments present in ALT slow .
We also compared differences in individual filament persistence during fast pacing. The total number of filaments during this pacing period was 527 and 3197 for ALT fast and ALT slow , respectively. Short-lived filaments dominated both ALT fast and ALT slow filament duration distributions (median duration was 17 ms for both, Fig. 5d), The duration of post-S2 arrhythmia was quantified by the number of times that reentry propagated around the pulmonary veins (referred to as a reentrant cycle). The period of reentries that sustained beyond one reentrant cycle was approximately 250 ms. When post-S2 reentrant activity lasted longer than three reentrant cycles, arrhythmia broke up into multiple wavelets (labeled "3+ "). The simulation indicated by the white cross had arrhythmia lasting longer than 2 s. and the difference in filament durations was not significant (Table 1). However, the variance in filament durations of the ALT fast versus ALT slow simulations was significantly different (1247.7 vs. 2432.6, AB = 525287.5, p = 0.003), with ALT slow exhibiting an increased proportion of long-duration filaments (Fig. 5d, top). Furthermore, 9 of the 14 long-lived filament outliers in the ALT slow simulation persisted longer than 400 ms, whereas only 1 of the 6 long-lived outliers did so in the ALT fast simulation. Thus, long lasting filaments occurred more frequently during fast pacing in atria with increased propensity to CDA. Regional differences in the density of points at which PSs occurred over time were observed during fast pacing. In ALT fast , there was higher phase singularity density in the LA than in the RA ( Supplementary Fig. S3). In ALT slow , phase singularity density was higher in the RA than in the LA except during 250-ms CL pacing, when the reverse was true. This was perhaps due to shorter APD in the LA (Fig. 1a), which allowed the LA to support more scroll waves.

S2 Stimulus Location
Arrhythmia persistence during fast pacing. New filaments can form when wavebreak creates new scroll waves, or when preexisting scroll waves interact or break up, resulting in the bifurcation and/or amalgamation of filaments. The proportion of filaments that formed from preexisting filaments during fast pacing (270-ms to 250-ms CL) was higher for ALT slow (0.57 vs. 0.51, χ 2 = 6.02, p = 0.014). Likewise, the proportion of filaments that were involved in forming new filaments was also higher for ALT slow (0.55 vs. 0.50, χ 2 = 3.87, p = 0.049). This suggested that increased filament interaction contributed to arrhythmia maintenance in ALT slow .
To further characterize the dynamics of filament formation and persistence, we analyzed filament trees (FTs), defined as the set of filaments that interact with each other spatiotemporally. During fast pacing, there were 198 and 976 total FTs in ALT fast and ALT slow , respectively. In both simulations, most FTs were short-lived (median duration was 37.5 ms and 52 ms for ALT fast and ALT slow , respectively) ( Fig. 5e), but FT duration tended to be longer in ALT slow ( Table 1). The distribution of FT durations in ALT slow was more right-skewed, with many more long-duration FT outliers (Fig. 5e, top).
Two examples of long-lived FT outliers from ALT fast and ALT slow are presented in Fig. 6a,b and Supplementary Videos S4 and S5. Although the temporal branching patterns of FTs could be quite complex, several FTs displayed an overall linear structure (Fig. 6a), indicating that each FT was dominated by a single scroll wave. The lifetimes of these scroll waves were often punctuated by wavebreak, resulting in the branching of a single filament into several short-lived filaments (< 25 ms). The ALT fast scroll wave meandered around the left atrial roof and lasted 540 ms (Fig. 6b, left). During its lifetime, one wavebreak occurred, separating the FT into two long-duration filaments (> 100 ms) occurring before and after the wavebreak (Fig. 6a, top). The ALT slow scroll wave meandered in a limited region on the right atrial wall near the superior vena cava and lasted 1,043 ms, substantially longer than the ALT fast FT (Fig. 6b, right). Four wavebreak episodes occurred, separating the ALT slow FT into five long-duration filaments  (> 75 ms, Fig. 6a, bottom). Thus, these patterns suggested that wavebreak contributed to increased persistence of scroll waves in ALT slow as compared to ALT fast . In order to evaluate how representative these FTs were of general differences in FT composition between ALT fast and ALT slow , we assessed four different FT composition metrics: mean filament duration, number of filaments, maximum filament duration, and sum of filament durations. The averages of each of these metrics over all FTs within a simulation, and statistical tests of the differences between simulations, are reported in Table 1. For all FT composition metrics, ALT slow tended to have larger values than ALT fast , with statistical significance in all metrics except mean filament duration. Not surprisingly, these metrics were also positively correlated with longer FT duration (Table 2). Figure 6c illustrates the relationship between FT duration and FT composition metrics. Although a modest linear trend was present in the relationship between FT duration and mean filament length, FTs with the longest duration deviated significantly from this trend, having very short mean filament duration (Fig. 6c, column 1). FT duration was more strongly correlated with the total number of filaments per FT (column 2), with the maximum filament duration within each FT (column 3), and with the sum of filament durations within each FT (column 4) (Supplementary Table S2).
These statistical trends were consistent with observations in the two example FTs from Fig. 6a,b. In general, the longest-lived FTs were composed of a few long-duration filaments (> 75 ms) preceded and/or followed by the numerous short-duration filaments (< 25 ms) that formed during wavebreak. Therefore, mean filament duration was weakly correlated with FT duration, but the number of filaments and the sum of filament durations were strongly correlated with FT duration (Table 2 and Fig. 6c). Additionally, maximum filament duration was more strongly correlated with FT duration in ALT fast than in ALT slow (Δ r = 0.13 [0.06, 0.17], Table 2). This was because long-lived FTs in ALT fast tended to be composed of only one or two long-duration filaments, while long-lived FTs in ALT slow tended to be composed of many long-duration filaments, due to the occurrence of several wavebreak episodes in each FTs lifetime. Thus, the above results demonstrate that increased wavebreak, resulting from increased propensity to CDA in ALT slow , helped to initiate and sustain arrhythmia in the remodeled atria during fast pacing.

Discussion
APD alternans has been observed at pacing rates as slow as 100-120 bpm in the atria of cardioverted AF patients undergoing electrophysiological study 6 . In a previous simulation study, we examined the underlying cellular mechanisms driving APD alternans in remodeled atrial myocytes, showing that altered RyR2 kinetics and down regulation of L-type calcium current induced by AF remodeling can lead to CDA at slower pacing rates 10 . Many remodeling changes secondary to AF are profibrillatory, leading to a cycle of worsening AF susceptibility over time (the idea that "AF begets AF"). In particular, shortened APD and decreased CV promote reentry and can stabilize scroll waves to help maintain AF 12 . Since alternans promotes breakup and thus destabilizes scroll waves, the persistence of the arrhythmias in the presence of alternans depends critically on restitution properties 13 . In this study, we sought to explore whether alternans occurring due to the calcium-dependent mechanisms we identified previously ultimately increases or decreases arrhythmia susceptibility.
Our results demonstrate that CDA occurring at slow pacing rates in remodeled atrial myocytes can lead to increased arrhythmia vulnerability, complexity, and maintenance. In a biophysically detailed, anatomically realistic model of the human atria, we found that shortened APD and reduced CV alone were not sufficient to render the atria vulnerable to reentry initiation by an ectopic beat from the pulmonary veins. In contrast, when CDA propensity was increased in the atria via the mechanisms we identified previously 10 , reentry could be induced at two of the right pulmonary vein locations tested. The atria model with increased CDA propensity also displayed more frequent wavebreak at fast pacing rates, resulting in more complex arrhythmia. Frequent wavebreak at fast pacing rates was also associated with longer-duration scroll waves, suggesting that increased wavebreak due to calcium handling instability ultimately increased arrhythmia persistence by producing more scroll waves that helped maintain the arrhythmia. These results suggest that elevated CDA propensity due to AF-induced remodeling of calcium handling may provide a basis for new antiarrhythmic strategies for AF. Three-dimensional computer models of atrial anatomy and electrophysiology provide a powerful tool for dissecting the mechanisms involved in AF arrhythmogenesis 14,15 . Previous computational studies have examined the arrhythmogenic role of decreased wavelength due to shortened APD and/or decreased CV [16][17][18] . Gong et al. 19 tested the effect of electrically driven alternans on AF initiation in a computer model of paroxysmal AF, in which upregulation of L-type calcium current caused in APD prolongation and steepening of APD restitution slope. However, APD and L-type calcium current are unchanged in paroxysmal AF patients 20 and are decreased in chronic AF patients 21 , so the clinical relevance of their findings may be limited. The cellular origins of APD alternans (e.g. steep APD restitution slope vs. calcium handling instability) are particularly important because the dynamical instabilities resulting from interaction of APD and CV restitution greatly affect spatial APD heterogeneity and hence arrhythmia vulnerability 11,22 . The mechanisms explored in this study stem from a more current understanding of the role of calcium handling abnormalities in human AF 10 and may thus provide better insight into the pathophysiological mechanisms underlying arrhythmogenesis in patients.
Our results demonstrate that in atria, increased propensity to CDA at slower pacing rates can significantly affect arrhythmia susceptibility. CDA created a vulnerable substrate defined by steep repolarization gradients, which allowed an ectopic beat to induce arrhythmia. This mechanism might explain clinical observations of premature beats initiating AF in patients during slow pacing rates when APD alternans was present 23 . Increased propensity to CDA at slower pacing rates also increased the complexity and persistence of arrhythmias occurring during fast pacing. Arrhythmias were characterized by unstable scroll waves which meandered over large areas and frequently broke up, similar to those observed clinically with ECG imaging 24 . Frequent wavebreak due to underlying calcium instability helped perpetuate arrhythmia in the remodeled atria, supporting the view that abnormal calcium handling may play a significant role in the AF substrate of patients with elevated alternans propensity.
The main mechanisms sustaining AF are much debated, with one of the central controversies being whether AF occurs primarily due to organized reentrant drivers or multiple wavelets 25 . Although several lines of evidence suggest that organized reentrant drivers underlie arrhythmias in many AF patients and may be targeted to effectively terminate arrhythmia 26 , patients with longer-duration continuous AF usually do not respond to ablation therapy 24 . Although exploring AF termination was beyond the scope of this study, we hypothesize that the disorganized, multiple-wavelet arrhythmia occurring due to increased CDA propensity in our model would be more difficult to terminate. Thus, CDA may reflect some of the progressive changes that occur in the underlying substrate during long-standing AF which make standard therapies ineffective. However, in the present study, we modeled fibrotic changes due to AF as a global decrease in tissue conductivity, without taking into account specific patterns of structural heterogeneity that can influence arrhythmia organization [27][28][29] . Further work is needed to determine how interaction of structural and functional mechanisms affect arrhythmia organization in AF 27,30,31 , particularly in patients that do not respond to ablation therapy.
Recent efforts to identify mechanisms and potential pharmacological targets in AF have shown that suppressing CDA may have significant antiarrhythmic benefits. Calcium transient alternans and abnormal calcium release events have been observed mouse 32 and canine 33 models of AF. Treatment with S107 in mouse 32 and dantrolene in canine hearts 34 appeared to stabilize RyR2 and effectively reverse abnormal changes in calcium handling to prevent arrhythmia. Our results suggest that similar mechanisms may be at play in human AF as well, since increased propensity to CDA driven by RyR2 dysfunction led to increased arrhythmia vulnerability, complexity, and persistence in our model. Recent evidence for this comes from experiments in human atrial myocytes of control and AF patients, which identified a link between adenosine signaling, alternans, and calcium handling instability, potentially mediated by PKA regulation of RyR2 35 . However, the mechanisms of RyR2 dysregulation in AF are extremely complex 36 and controversial 37 . Further work needs to be done to clarify the mechanisms of RyR2 regulation and to explore the antiarrhythmic benefits of RyR2-targeting therapies that address the dual role of calcium as both an initiating trigger and vulnerable substrate in AF 5 .
Although the present study focused on the role of alternans in chronic AF, elevated propensity to APD alternans has been observed in paroxysmal AF patients as well 23,38,39 . However, the mechanisms underlying alternans may differ, since calcium handling is altered between paroxysmal and chronic AF 20 . In paroxysmal AF, evidence suggests that SK channels 40 , acetylcholine-activated K + current 41 , and late Na + current 42 may play a more prominent role in alternans propensity and arrhythmogenesis. Given the arrhythmogenic implications of APD alternans presented in this and other studies, these differences in the fundamental mechanisms underlying APD alternans in paroxysmal versus chronic AF underscore the need for further research on the diversity of AF mechanisms 2 , in order to develop better targeted, patient-specific therapies.

Methods
Human atria model. We used an anatomically realistic computer model of the human atria to investigate the effect of increased CDA propensity on arrhythmia initiation and maintenance in the atria. Model geometry was based on the Visible Human female dataset 43 , which has been used previously to investigate mechanisms of arrhythmogenesis related to AF 17 . To explore the effects of CDA secondary to AF-induced electrical remodeling on arrhythmogenesis, we tested two atria with chronic AF-induced structural and electrical remodeling. Structural changes were modeled by scaling anisotropic heterogeneous tissue conductivities to produce an average CV of 40 cm/s, consistent with slowed conduction in chronic AF 44 (Supplementary Table S1). Electrical remodeling was based on a human atrial AP model described by Grandi et al., which produced shortened APD as observed in chronic AF 45 . These parameters led to alternans onset at fast pacing rates in the first atria model, ALT fast . The second atria model, ALT slow , incorporated an additional remodeling feature: reduction of RyR2 inactivation rate by 50%, which we have shown previously results in increased propensity to CDA, matching alternans onset at slower pacing rates observed clinically 10 (Fig. 1a, see Supplementary Methods). Data on ionic heterogeneity in the atria during AF is limited, so only regional AP heterogeneity between the left (Fig. 1b) and right (Fig. 1c) atrium were included in the atria models (Fig. 1d), as described by Grandi et al. 45 .
Numerical methods. A monodomain formulation was used to represent electrical propagation throughout the human atria model 46

Arrhythmia induction protocols.
To compare arrhythmia dynamics between ALT fast and ALT slow models, we performed a dynamic pacing protocol similar to one used clinically to assess the presence of alternans and induce AF in patients 23 . All nodes in the model were initialized using steady-state values from a single cell paced at 750-ms CL. The atria were paced from the SAN region at 750-ms CL for 20 beats (Fig. 1e), and then for 32 beats at each subsequent pacing CL: 500-ms to 300-ms CL in 50-ms decrements, then 290-ms to 250-ms CL in 10-ms decrements. Arrhythmias were examined (see Methods section on filament dynamics) during the fastest pacing rates when wavebreak and reentry occurred (CL ≤ 270 ms).
To assess arrhythmia vulnerability due to ectopic beats at slower pacing rates, we performed an S1-S2 pacing protocol using the same initial dynamic pacing protocol described above. After the 16 th even or odd beat at 300-ms CL pacing, SAN pacing ceased and an ectopic stimulus (S2) was applied to one of eight locations in pulmonary veins at eleven different CIs, ranging from 250 ms to 300 ms. Figure 1e shows location of the S2 stimuli, in the right superior pulmonary veins (RS0-3) or the right inferior pulmonary veins (RI0-3). Post-S2 activity was simulated for 2 s or until activity terminated, whichever came first.
Quantification of APD alternans. APD was calculated as the time from maximal upstroke velocity to 90% repolarization of the transmembrane potential (V m ) from phase II amplitude. Alternans magnitude was quantified as the mean change in APD from beat to beat over the last 10 pairs of beats (11 beats total) at a given pacing CL 23 .

Filament dynamics.
Phase singularities are points where the phase of the AP is undefined, occurring at locations of wavebreak or at the tip of a reentrant spiral wave. In 3 D, phase singularities form the vortex filaments of reentrant scroll waves. To robustly detect filaments in our simulations despite variations in AP morphology 49 , we developed a novel method for filament localization, which involved constructing phase maps from activation times (see Supplementary Methods). To gain insight into filament dynamics, we tracked at 1-kHz resolution each filament's birth, death, and interactions with other filaments, splitting (bifurcation) or merging (amalgamation) to form new filaments. A FT was defined as the set of filaments that interact with each other over time. Spurious filaments associated with brief wavebreak were not included in the dataset to allow interpretation of the temporal branching structure of filament interactions. An example snapshot of filament localization during arrhythmia is shown in Fig. 1f. Details on filament detection can be found in Supplementary Methods. Statistical analysis. Adjusted box-whisker plots were used to compare skewed distributions and detect outliers 50 . Chi-squared tests of independence without correction were used to assess differences in filament interaction 51 . The Mann-Whitney U statistic was used to test for differences in distribution location. The Ansari-Bradley statistic (AB) was used to test for differences in the scale (i.e. variance) of distributions with identical medians 52 . Two-tailed tests were performed with p < 0.05 considered significant. Pearson's correlation coefficient r was used to quantify the linear dependence between variables. Significance of Pearson correlations and their differences (Δ r) was assessed using a percentile bootstrap method to estimate 95% confidence intervals (Wilcox-Muska confidence intervals) in the independent case, and using the HC4 method in the overlapping dependent case 53 . These methods have been shown to perform well under conditions of nonnormality and heteroscedasticity 54 .