The human whole-atria anatomical model described in Seemann et al. (2006) was used as the basis to construct the six virtual human patient models of AF. The whole-atria anatomical model includes fiber orientation, conduction anisotropy, and spatial heterogeneities in ionic currents and conduction velocity (CV) in the main atrial structures: LA, RA, sinoatrial node (SAN), cresta terminalis, pectinate muscles, fossa ovalis, Bachmann’s bundle, cavotricuspid isthmus, left atrial appendage (LAPG), right atrial appendage (RAPG), atrioventricular ring, and interatrial bridges (Seemann et al., 2006; Krueger et al., 2011; Tobon et al., 2013). Anisotropic ratio (transversal to longitudinal ratio of conductivity) and spatial heterogeneities are summarized in Table S1 in Supplementary Material.

Each human whole-atria model was assigned a different AP phenotype with human atrial membrane kinetics based on the Maleckar et al. model (Maleckar et al., 2009), but with ionic conductances calibrated to account for AF remodeling and intersubject variability in APD using an experimentally calibrated population of human AP models as described in Sánchez et al. (2014). In short, a population of 2,275 AP atrial models was constructed based on an AF-remodeled version of the Maleckar et al. model (70% decrease in I_CaL, 50% decrease in I_to, 50% decrease in I_Kur, and 100% increase in I_K1), considering variability in six ionic current conductances (I_K1, I_NaK, I_CaL, I_to, I_Kur, and I_NaCa) and calibrated using AP recordings from human atrial trabeculae (Wettwer et al., 2004, 2013). In order to evaluate effects of differences in APD₉₀, APD₅₀, or APD₂₀ values on atrial dynamics, the human AP models (Figure 1A) were classified into six subpopulations based on whether they exhibited long or short APD₉₀, APD₅₀, or APD₂₀ values (Figures 1B–D). For each degree of repolarization (APD₉₀, APD₅₀, or APD₂₀), the human atrial models yielding APD within the first quartile in the overall population were classified as “short APD,” and the models above the third quartile as “long APD.” For example, the short APD₉₀ subpopulation included all models with APD₉₀ shorter than 172.6 ms and the long APD₉₀ subpopulation those with APD₉₀ longer than 222.3 ms (first and third APD₉₀ quartiles, respectively, Figure 1B).

In human atrial cardiomyocytes, APD₉₀, APD₅₀, and APD₂₀ are frequently (although not always) related with each other, i.e., long APD₉₀ is in general associated with long APD₅₀ and long APD₂₀, as shown in Figure 1E via scatter plots of experimentally measured APD₉₀ vs APD₅₀, APD₉₀ vs APD₂₀, and APD₅₀ vs APD₂₀. The scatter plots of Figure 1E further illustrate natural overlaps between short/long APD categories (e.g., long APD₉₀ vs short APD₅₀; Figure 1E, left panel) and less likely intersections (e.g., long APD₂₀ vs short APD₅₀; Figure 1E, right panel). Such overlaps are also captured in the six generated model subpopulations as illustrated in the Venn diagrams displayed in Figure 1F, showing frequency of intersections in number of APs between the different (and more largely sampled) APD subpopulations.

The distributions of ionic current conductances in each of the short/long APD₉₀, APD₅₀, and APD₂₀ subpopulations are shown in Figure 2. This is in agreement with experimental studies showing variability in cellular electrophysiology in human cardiomyocytes and, specifically, in the atria (Wang et al., 1993). Each subpopulation of human AP models was used to assign membrane kinetics in one of the six human whole-atria models, by randomly assigning an AP model of the corresponding subpopulation to each cell (once per whole-atrial model). Figure S1 in Supplementary Material shows APD distributions in the six whole-atria models in sinus rhythm. These distributions indicate that the effect of APD heterogeneity due to the existence of different atrial structures has a stronger impact on APD dispersion than the random APD variability within each subpopulation. This means that tissue coupling and macrostructural regional differences in the anatomical model (Table S1 in Supplementary Material) attenuate cell-to-cell differences in APD due to random sampling from the AP subpopulations.

Our simulations also incorporate the effect of a concentration of 1 nM acetylcholine ([ACh]) in the whole models, simulated as in previous studies (Kneller et al., 2002), in order to facilitate the inducibility of arrhythmic behavior (Pandit et al., 2011) and to incorporate electrophysiological heterogeneities in the different atrial structures (Table S1 in Supplementary Material). Given the uncertainty on exact distributions of ACh release in the human atria (Kneller et al., 2002; Jones et al., 2012), we adopted a simplified approach for its modeling, accounting for left and right atrial (LA/RA) differences on I_KACh as reported in Sarmast et al. (2003). Larger ACh concentrations yielded an overall reduction in the regularity indices described below, with similar spatial distributions (Figure S2 in Supplementary Material).

A preliminary conditioning of each whole-atria model was performed by applying a train of 40 periodic stimuli at the SAN site at a cycle length of 500 ms. Then, six extra-stimuli were applied near the right pulmonary veins (RPV) in the LA to generate reentries and fibrillatory behavior. The cycle length of these six extra-stimuli was different for each AP phenotype (in the range 130–180 ms). All externally applied stimuli were 2-ms duration and twice diastolic threshold current pulses.

In order to investigate how APD differences may modulate the response to antiarrhythmic therapies, the effects of ionic current block on fibrillatory dynamics and reentrant properties were quantified. Given their importance in atrial dynamics, as identified in previous studies (Kneller et al., 2005; Pandit et al., 2005, 2011; Sánchez et al., 2012), 30% I_K1 block, 30% I_NaK block, and 15% I_Na block were simulated in all six whole-atrial models. These degrees of block allow the analysis of fibrillatory properties without critical alterations in cellular electrophysiology and/or electrical propagation (Sánchez et al., 2012), in agreement with experimental studies on the inhibition of I_K1 with barium (Wu et al., 1999), I_Na with flecainide (Barekatain and Razavi, 2012), and I_NaK with digitalis (Wasserstrom and Aistrup, 2005). Furthermore, similar degrees of block have been widely used in sensitivity analysis studies (Romero et al., 2009, 2011; Keller et al., 2010; Sánchez et al., 2012; Chang et al., 2015). Simulation of these ionic blocks was initiated 3 s after arrhythmia generation in each of the six AF virtual models.

The EGM indices obtained from the six human whole-atria models were compared to intracardiac recordings from AF patients included in the Ann Arbor database (Ann Arbor EGM Libraries, Chicago, IL, USA) (see Figures S3 and S4 in Supplementary Material). We considered only the unipolar EGMs registered on the atrial surface, i.e., those in the high RA and the RAPG. The sampling rate for signal acquisition was 1,000 Hz. The signal processing techniques described in the following section and the Supplementary Material for the simulated EGMs were also applied to calculate dominant frequency (DF), organization index (OI), and regularity index (RI) in the EGM database. Further details in metrics calculation are provided in the Supplementary Material. Coupling index (CP) was not calculated as only one EGM signal per patient was available (Faes and Ravelli, 2007).

In the six virtual models, EGMs were obtained in 49 electrodes covering both the LA and the RA and located at a distance of about 2 mm from the endocardial surface, as in unipolar EGM recordings (Plonsey and Rudy, 1980; Gima and Rudy, 2002; Baher et al., 2007). A fourth-order Butterworth band-pass-filter (f_c1 = 40 Hz; f_c2 = 250 Hz), followed by a signal-rectifier and a fourth-order Butterworth high-pass-filter (f_c = 20 Hz), were applied to the EGMs in order to minimize spectral components different from those of the main activation wavefront (Botteron and Smith, 1995). The four indices mentioned in the previous paragraph, DF, OI, RI, and CP, were extracted from the 49 EGM signals and interpolated to the whole-atrial tissue for representation.

Whole-atrial simulations were conducted using the ELVIRA code, which solves the monodomain equation using the finite element method as described in Heidenreich et al. (2010), with a time step of 0.04 ms ensuring numerical convergence of the results (Figure S5 in Supplementary Material) and minimizing computation time. Signal processing techniques were implemented using Matlab.

Figure 3 illustrates the strong effect of differences in APD in the dynamics of human AF through quantification of DF, RI, and OI in the six virtual patients. AF dynamics were modulated by differences in all three phases of repolarization, as shown in Figure 3 for the models with short vs long APD₉₀, APD₅₀, and APD₂₀. Furthermore, the computed EGMs in the six whole-atria models yielded all arrhythmia-related indices in agreement with those in the Ann Arbor EGM database, as shown in the Supplementary Material.

Fibrillatory dynamics in the human long APD₉₀ whole-atria model was sustained by a single stable rotor in the RA, with slower DF than in the short APD₉₀ atrial model (top-left panel in Figure 3). The short APD₉₀ model exhibited more chaotic AF dynamics (sustained by several reentrant rotors) with high DF values, particularly in the LA and low RI (reflecting changes in activation patterns in different cycles). OI was generally high for both phenotypes, meaning that reentrant rotors presented no frequential components significantly different from DF and its harmonics.

The human model with long APD₅₀ exhibited smaller values of DF (1.77 Hz lower in average) and OI (0.42 smaller in average) in the RA than the model with short APD₅₀ (second row-left panels in Figure 3). The main reason was a block of interatrial propagation occurring in one out of three wavefronts from the LA to the RA for the long APD₅₀ simulation. However, RI and CP were similar for both human models, suggesting local morphological changes of similar magnitude in the EGMs over time and between neighboring regions (second row-right panels in Figure 3).

Interestingly, AF dynamics in human atrial models with short vs long APD₂₀ displayed notable differences in fibrillatory patterns. The human atrial long APD₂₀ model was associated with very organized and regular activation patterns (OI and RI indices) as compared with those in the human atrial short APD₂₀ model. AF dynamics in the short APD₂₀ model were similar to those in the short APD₅₀ model but with a more chaotic propagation in the RA due to the wavefront meandering, as shown by the lower values of both OI and RI (second and third rows, middle panels in Figure 3). Regarding the long APD₂₀, there was only one main rotor, with low DF, high OI, high RI, and high CP in all the tissue (third row in Figure 3), which stopped propagating after 6.3 s of simulation due to interatrial conduction block.

In spite of differences in atrial dynamics shown in Figure 3, ionic current block yielded similar propagation patterns in all six virtual whole-atria models, as illustrated in Figure 4. The cores of the main reentrant rotors were usually located near the RPV, the LPV, the LAPG, and the RAPG, and secondary rotors usually appeared close to the SVC and the MV. Some particular cases, such as I_K1 block in the short APD₅₀ virtual phenotype as well as I_NaK block in both short and long APD₅₀ virtual phenotypes, presented significant wave meandering and temporary rotors in the RA. In 17 out of the 18 considered interventions, the fibrillatory behavior remained self-sustained for the duration of the simulation. Only 30% I_NaK block in the long APD₉₀ virtual model led to arrhythmia termination. Complete information about the number of rotors and reentrant circuits, i.e., those without phase singularities but rotating around an anatomical structure, in the LA and the RA are displayed in Table S2 in Supplementary Material classified by duration and path length for each simulation. The antiarrhythmic effect of the ionic current blocks simulated was, therefore, limited in the simulations, even though the complexity of the arrhythmia was decreased.

Overall, I_K1, I_NaK, or I_Na block led to smaller DF and more uniform DF distributions in the six virtual phenotypes (Figure 4). This was due to the predominance of the main rotors over secondary rotors entailing additional periodicity of atrial tissue activation. However, significant LA–RA gradients in DF arose in some scenarios due to interatrial propagation block. This was the case in simulations with I_NaK and I_Na blocks in the short APD₂₀ virtual phenotype (third row, right panels in Figure 4). Both I_K1 block in the short APD₅₀ phenotype and I_NaK block in the long APD₅₀ phenotype entailed low DF in most atrial tissue except for high DF in reduced areas close to the junctions between the LA and the RA (second row, middle panels in Figure 4). This regional high DF was produced by a different delay in the wavefronts propagating and reentering in each atrium due to the APD alterations derived from the ionic blocks. Two depolarizations per rotor cycle were observed in the two cases: one due to the propagation wavefront rotating in the RA and the other due to the income of an earlier wavefront from the LA, colliding afterward with the RA wavefront. These areas did not have time to repolarize completely after their activation by each wavefront and, therefore, depolarizations were locally less prominent.

Regarding OI distributions, OI was higher in the LA than in the RA in all the simulated conditions due to the wave meandering and generation of secondary wavelets in the RA. This, thus, led to more disorganized activation patterns in RA than in LA (Figure 5). A significant reduction of OI in the LA was only observed after inhibition of I_K1 in the short APD₂₀ phenotype (third row in Figure 5). The highest OI values were obtained when I_K1, I_NaK, or I_Na were inhibited in the long APD₉₀ phenotype, although these were slightly smaller than those in the phenotype when no ionic blocks were simulated (first row in Figure 5).

Figures 6 and 7 show results of RI distributions and CP between adjacent electrodes, respectively. RI values after ionic inhibition were all smaller than in their respective default cases. Furthermore, they were generally higher in the LA than in the RA, highlighting important LA–RA gradients in most of the simulated cases because of the presence of RA rotors with acute meandering, whereas LA activation was usually more stable. As for OI, inhibition of I_K1 in the short APD₂₀ phenotype led to very irregular activation patterns, i.e., low RI, in both atria (third row, second column in Figure 6). The most regular activation patterns were obtained following I_Na block in the long APD₉₀ phenotype (top-right panel in Figure 6). Regarding CP, a decrease in coupling with respect to the default cases was generally observed after the simulations of ionic blocks (Figure 7). The strongest couplings were observed in the LA after inhibiting I_Na in four of the six virtual phenotypes, whereas they were very weak in both atria for the short APD₅₀ and APD₂₀ phenotypes after blocks of any of the three ionic currents (Figure 7).

Percentual variations of the analyzed indices following the three ionic blocks with respect to the averaged values in the six virtual phenotypes are shown in Figure 8 for both the LA (top panels) and the RA (bottom panels). Bars of the same color in Figure 8 represent each virtual phenotype, whereas the analyzed indices are shown in the horizontal axis. The effects of these ionic blocks were qualitatively similar between virtual phenotypes, except for some particular cases. For example, the block of any of the three ionic currents promoted interatrial differences in RI and CP in the short APD₉₀ phenotype by increasing both indices in the LA and decreasing them in the RA. The opposite occurred in the short APD₂₀ phenotype, in which RI and CP were decreased in the LA and increased in the RA.

In this study, the Maleckar et al. model (Maleckar et al., 2009) was used to simulate cellular electrophysiology, given its improved representation of human atrial repolarization dynamics. The choice of a different human atrial cell model could make results be different from those shown in this study. Nevertheless, AP phenotypes in AF with other human atrial cell models (Courtemanche et al., 1998; Grandi et al., 2011) after experimental calibration were shown to be similar to those with the Maleckar et al. model in a previous study (Sánchez et al., 2014).

Reentrant behavior was generated using a particular stimulation protocol consisting of periodic stimuli applied at the SAN followed by six extra-stimuli near the RPV. However, AF in clinical practice may present multiple ectopic foci that trigger and sustain the arrhythmia, and their location may be variable as well (Jalife et al., 2002). Furthermore, a healthy atrial anatomical model was used in this study, in order to focus on the effects of AF-induced electrophysiological rather than structural remodeling (like fibrosis), which would also affect arrhythmia generation and maintenance (Burstein and Nattel, 2008).

The effects of selective ionic blocks were analyzed individually in this study. However, many antiarrhythmic drugs, such as flecainide, digitalis, or barium, have multiple ionic targets and their simultaneous inhibitions could have synergistic and/or non-linear effects on arrhythmic markers. Therefore, a multi-target approach like that used in Sivagangabalan et al. (2014) could be interesting to explore in future studies. In addition, we only simulated moderate selective ion channel inhibitions (30% I_K1 block, 30% I_NaK block, and 15% I_Na block). Larger ion channel blocks might lead to a successful termination of atrial fibrillatory dynamics in more virtual patients, but they would also have stronger effects in the ventricles, raising potential safety concerns.

The relatively small number of virtual electrodes used (49) led to visible interpolation artifacts in the figures. Nevertheless, this precision was high enough to allow proper characterization of arrhythmic behavior, since the electrical propagation wavefronts were much wider than the separation between adjacent electrodes.

As a final remark, the random assignment of a specific model within each subpopulation to each cell in the 3D model could potentially lead to neighboring cells with very different electrophysiological properties. This effect was, however, reduced for two reasons: (1) variability in the AP within each subpopulation was relatively small (first and fourth quartiles of the subpopulations were chosen to represent the slight variability present within one subject); (2) the electrotonic coupling between cells strongly smoothed electrical gradients during propagation (Figure S1 in Supplementary Material).

Inter-patient variability in cell repolarization entails notable differences in the organization and stability of arrhythmias in AF-electrophysiologically-remodeled tissue. Specific ionic blocks of I_K1, I_NaK, and I_Na reduce the regularity of fibrillation patterns and promote interatrial differences, especially for the phenotypes with short APD at 20, 50, and/or 90% repolarization.