The Major Role of I K1 in Mechanisms of Rotor Drift in the Atria: A Computational Study

Abstract

Maintenance of paroxysmal atrial fibrillation (AF) by fast rotors in the left atrium (LA) or at the pulmonary veins (PVs) is not fully understood. This review describes the role of the heterogeneous distribution of transmembrane currents in the PVs and LA junction (PV-LAJ) in the localization of rotors in the PVs. Experimentally observed heterogeneities in I_K1, I_Ks, I_Kr, I_to, and I_CaL in the PV-LAJ were incorporated into models of human atrial kinetics to simulate various conditions and investigate rotor drifting mechanisms. Spatial gradients in the currents resulted in shorter action potential duration, less negative minimum diastolic potential, slower upstroke and conduction velocity for rotors in the PV region than in the LA. Rotors under such conditions drifted toward the PV and stabilized at the less excitable region. Our simulations suggest that I_K1 heterogeneity is dominant in determining the drift direction through its impact on the excitability gradient. These results provide a novel framework for understanding the complex dynamics of rotors in AF.

Introduction

The mechanisms of atrial fibrillation (AF), the most common cardiac arrhythmia in the clinical practice, are not fully understood. Acute AF in normal isolated sheep hearts has been found to often depend on fast rotors localized mainly to the posterior wall of the left atrium (LA) and the pulmonary veins (PVs) junction (PV-LAJ) with fibrillatory conduction toward the rest of the atria. ¹ Recent clinical data also point to rotors in various atrial sites as a mechanism driving paroxysmal AF.^2–4 Previous simulations have already demonstrated that rotors can occur in the PV-LAJ, provided the PV size is adequate, and nonuniform coupling conditions exist. ⁵ However, how the rotors form or drift in the PV-LAJ, as well as the drift underlying ionic mechanisms, have not been investigated systematically.^6,7 Earlier studies using cardiac computer models have shown a causal link between rotor drift and spatial heterogeneity in the action potential (AP) properties. ⁸ Rotors drifted toward regions with prolonged action potential duration (APD) or reduced excitability,^9,10 which were inscribed mainly due to heterogeneities in multiple ion channels, ¹¹ and particularly K⁺ channels. ¹⁰

As shown in Figure 1, I_K1 and I_Kr play a role in rotor dynamics by affecting both the membrane APD and excitability.^12,13 Figure 1A presents the average duration and frequency of the tachyarrhythmic episodes in the wild-type (WT) and transgenic (TG) mice, in which I_K1 was upregulated by overexpressing the Kir2.1. ¹³ There was statistical difference between duration (WT, 3 ± 9; TG, 350 ± 1181 seconds, P < 0.001) and frequency (WT, 26.2 ± 5.2; TG, 44.6 ± 4.3 Hz, P < 0.001) of the arrhythmias in the two groups. In addition, Figure 1A displays a representative snapshot of the distribution of membrane voltage and the underlying inactivation gating variables (h·j) of the sodium current (I_Na) during simulated reentry in WT and TG mouse hearts. The product of the h and j gates is representative of the excitability of the resting tissue. The voltage maps show that the wavelength is shorter and the diastolic potential is more hyperpolarized (darker blue) over a wider area in the TG cardiac sheet, compared to the WT (values of diastolic potentials were ~-94 mV in TG versus ~-90 mV in WT). This results in greater recovery of Na⁺ channel availability (maximum value of h·j is 0.808 versus 0.684 in TG versus WT hearts, respectively). Furthermore, the accelerated final phase of repolarization results in a larger area in which Na+ channels have recovered (orange), and thus greater excitability in front of the depolarization wave front. These factors contribute in part to the higher rotor frequency seen in the TG hearts.

Figure 1.

The effect of I_K1 and I_Kr on rotor stability and frequency. (A) Top: average duration and frequency of arrhythmias in four WT and eight TG hearts from mice. Bottom: spatial distribution of transmembrane voltage (top) and inactivation gating variables (h·j gates, bottom) of the I_Na during stable reentry in a 2 × 2 cm² cardiac cell sheet. The diastolic potential is more hyperpolarized over a wider area in the TG heart sheet (~-94 mV), compared to the WT case (~-90 mV). The h·j gate map shows that the maximum availability of Na+ channels in TG (0.808) exceeds that of the WT heart (0.684) and occupies a larger area during the excitable gap (from Ref. 13). (B) The effects of hERG infection on reentry frequency in neonatal rat ventricular myocytes monolayers. Top: snapshots from representative optical mapping movies in monolayers infected with Ad-GFP (left), Ad-hERG (middle), or Ad-G628S (right). Bottom: reentry frequencies in individual monolayers infected with Ad-GFP (I_Kr-overexpression, open circles), Ad-hERG (P < 0.05 versus GFP, closed circles), or Ad-G628S (P = NS versus Ad-GFP; P < 0.05 versus Ad-hERG, open triangles). Horizontal bars represent mean values. One-way analysis of variance (ANOVA) with Tukey's multiple comparison test (reproduced from Ref. 12).

In addition, sustained functional reentry is seen to accelerate in cardiomyocytes' monolayers infected by hERG to increase their I_Kr. ¹² The top of Figure 1B shows phase maps of single rotors generating spiral waves in control (Ad-GFP, left), I_Kr-overexpressing (middle), and G628S-overexpressing (I_Kr blockade, right) monolayers of neonatal rat ventricular myocytes. ¹² The bottom graph in Figure 1B compares the ranges of individual rotation frequencies of control monolayers with that of I_Kr-overexpressing and G628S-overexpressing monolayers. The important role of I_Kr is illustrated by the significantly higher average rotational frequency in I_Kr-overexpressing monolayers than in GFP controls (21.12 ± 0.81, n = 43 versus 9.21 ± 0.58 Hz, n = 16; P < 0.001), whereas the frequency in the G628S monolayers is slightly lower (6.14 ± 0.3 Hz, n = 17; P = not significant [NS] versus GFP).

While the effect of the individual I_K1 and I_Kr on rotor dynamics and hence on its drift should be discernable, the complex coexisting heterogeneity in the expression of channels at the PV-LAJ, such as in dogs, which show a larger current density of I_Kr and a smaller density of I_K1 in the PV compared to LA,^14,15 in conjunction with structural discontinuities (narrow PV sleeves compared to the larger LA), precludes a simple prediction of the rotor dynamics at the PV-LAJ. To address such substrate effects, the data presented in this paper aim to elucidate ionic mechanisms of rotor drifting at the PV-LAJ, a region that is important for paroxysmal AF maintenance in many patients. We review here an investigation in two-dimensional (2D) and pseudo-three-dimensional (3D) models of the PV-LAJ to test the hypothesis that the characteristic heterogeneous dispersion of transmembrane currents during paroxysmal AF, in the absence of remodeling, is a determinant of rotor drift. Our numerical simulations further demonstrate I_K1 dominance in conveying a preferential rotor drift direction toward the PVs. We also demonstrate the cycle-by-cycle mechanism by which PV-LAJ regions with longer refractoriness and lower excitability tend to attract rotors. Notwithstanding the reliance on simple models and animal-based ionic kinetics, it is proposed that the specific ionic properties of the human atrial substrate also predispose the PVs, or any other atrial region with reduced I_K1 density, to attract rotor activity during AF.^16–18

Numerical Methods

Computer simulations were performed on two simple models of the junction between the LA and the PV (Fig. 2) ¹⁹ :

1. 1.

   A 50 × 50 mm² regular 2D square mesh model was implemented and subjected to no-flux boundary conditions.

2. 2.

   A pseudo-3D cylindrical surface model of 50 mm length and perimeter constructed by applying no-flux boundary conditions at the LA and PV edges and periodic boundary conditions on the other two edges of a regular 2D mesh.

Figure 2.

The PV-LAJ models. (A) Conductances for I_CaL, I_Kr, I_Ks, I_to, and I_K1 follow a Boltzmann function (X_1/2 = 25 mm, vertical dashed line; ΔX = 5 mm) from LA to PV. Color bar: normalized Boltzmann factor across the junction. (B) Simulations were performed on (1) 2D model of 50 × 50 mm² and (2) cylindrical pseudo-3D model generated by a periodic boundary condition in the 2D model. Cyan-yellow is the steepest gradient region in the PV-LAJ (reproduced from Ref. 19).

The transmembrane potential was based on the Courtemanche-Ramirez-Nattel and Kneller (CRN-K) model of human atrial cell kinetics for AF without remodeling and in the presence of 0.0015 μM acetylcholine.^20–22 Electrical activity was computed using a monodomain and numerically stable Euler forward scheme with 100 μm internodal distance, 0.005 ms time step, and isotropic diffusion coefficient of 0.062 mm²/ms adjusted for conduction velocity of 48 cm/s.^23,24 Initial conditions were the membrane variables at resting potential following pacing a single cell at 1 Hz for 10 seconds.

Heterogeneous ionic conditions were implemented by assigning a spatial Boltzmann distribution of conductance for I_K1, I_Ks, I_Kr, I_to, and I_CaL between the LA and the PV based on the data from dogs ¹⁴ and described in Figure 2A (maximum conductance for I_Kr and I_Ks was increased by 50% and 60%, respectively, while that for I_K1, I_to, and I_CaL was diminished by 42%, 25%, and 30%, respectively). Single and paired ionic currents heterogeneities were implemented, as well as three additional conditions of ionic heterogeneity were modeled: (1) all currents varied spatially according to their corresponding Boltzmann function (Condition I); (2) all currents, except I_K1, varied (Condition II); and (3) only I_K1 varied in space (Condition III). Reentrant excitation patterns were induced by S1-S2 cross-field stimulation protocols, ²² and their pivoting points were identified as a singularity point (SP) and tracked in the phase domain from the moment of their appearance. ¹⁹

Spatial profiles of peak sodium channels availability (h·j)_peak were obtained by averaging those parameters in time for every pixel across the model in the last 2 of 5 seconds simulations. Measurements of minimum diastolic potential (MDP) and maximum upstroke velocity (dV/dt_max) followed similar procedure.

Results

Ionic Gradients and Rotor Attraction toward the PV

We tested the effect of all individual and pairwise currents on drift by switching them from heterogeneous to uniform (Fig. 3). Figure 3A shows that different ionic current heterogeneity, when applied individually, has a different drift direction. However, when pairs of ionic heterogeneities were combined, Figure 3B shows that the I_K1 heterogeneity is the primary ionic factor that determines the direction of the drift toward the PV; whatever I_K1 heterogeneity is combined with, the drift direction is always in the direction imposed by the I_K1, namely a drift from the LA toward the PV region.

Figure 3.

Spatial gradients of individual and paired currents and reentry drift in a flat model. (A) Effect of spatial gradient for each individual current on rotor drift. 2D PV-LAJ model showing spatiotemporal trajectories of SPs of rotors in gradients of I_Ks (black), I_Kr (green), I_to (yellow), I_CaL (blue), and I_K1 (red) as described by the individual Boltzmann distributions in Figure 2A. While I_to and I_K1 cause the rotor to drift toward the PV edge, the gradients in the other currents cause the rotor to drift toward the opposite LA edge. Insert: top view of the trajectories. (B) Effect of spatial gradients in paired currents on rotor drift. The table shows the direction of the drift when ionic gradients as indicated in the abscissa and ordinate (green cells) were considered. I_K1 is the only current whose gradient leads to PV attraction (in red fonts) of the rotor, regardless of gradient in any other current (reproduced from Ref. 19).

As a demonstration of major dependence of rotor drift in the PV-LAJ on the specific ionic current heterogeneity of I_K1 in contrast to all the other combined heterogeneities, we generated rotors in the three ionic conditions implemented in the cylindrical and flat model and tracked the spontaneous trajectory of their SP. Figure 4 demonstrates that the rotor in the heterogeneous PV-LAJ region drifts predictably as a consequence of the particular condition simulated.

Figure 4.

Simulations of rotor drift in the PV-LAJ and I_K1 role against all other currents. (A) Rotor dynamics on the cylindrical model using the CRN-K kinetics under three conditions: (1) all currents varied spatially, (2) all currents varied spatially, except I_K1, and (3) only I_K1 varied spatially. Color-coded traces show the SP trajectory of rotors initiated at the middle of the models (blue dots). Rotors under Conditions I and III drifted toward the PV edge. Drift direction was reversed when I_K1 heterogeneity was excluded in Condition II. (B) Time–space plots and arrows show rotor drift for the three conditions in a flat 2D model (reproduced from Ref. 19).

In Condition I, ie, heterogeneity in the currents I_K1, I_Ks, I_Kr, I_to, and I_CaL as characterized at the PV-LAJ in dog, ¹⁴ the drift is toward the PV edge of the model. However, when all the currents except I_K1 are set to be heterogeneous, and I_K1 density is maintained homogenously and is equal to the LA region value (Condition II), the drift reverses toward the LA (when the I_K1 is homogeneously equal to the PV region – not shown – the drift is also toward the LA). On the other hand, when all the currents are set uniform (either with the LA or the PV values) and only I_K1 is set to disperse as in the dog (Condition III), the drift is again toward the PV, but with a faster rate as compared with the drift in Condition I. This set of three scenarios in cylindrical as well as in the flat models, shown in Figure 2, clearly points to the strong effect of I_K1 on the direction of the rotor drift in the LA-PV junction area. As controls, we simulated rotors in 2D models with uniform ionic properties of either the LA or the PV that revealed nondrifting rotors. The rotors in the model with uniform LA properties were slightly faster than those in the model with the uniform PV properties (7.7 versus 7 Hz, respectively; other properties of the rotors did not vary by more than about 10%). ¹⁹

Heterogeneous Excitability and Rotor Drift in the PV-LAJ

Ionic heterogeneity is imposing nonuniform excitability properties and a rotor drift toward a predictable direction. ⁹ In Figure 5, we analyze the role of tissue excitability in rotor drifts by quantifying the spatiotemporal distribution of the product of membrane model parameters h and j (h·j), which determines the I_Na availability during rotor pivoting. ²² Figure 5A illustrates the drift trajectory of a rotor in a 2D model under Condition I and shows a time–space plot (TSP) of h·j along the line of the drift. A closer look at the TSP near the pivoting location during a single cycle shows higher values of h·j at the left side (LA) of the drift than at the right side (PV). Thus, a progressive shift of the pivoting point toward the PV is always toward a region with lower h·j values.

Figure 5.

Heterogeneous sodium availability during rotor drift in a flat model. (A) Left: voltage snapshot (gray scale) with superimposed trajectory of rotor SP (yellow tracing) drifting from the center of the model (red dot) toward its final location at the PV edge (blue dot). Right: color coded time–space plot for h·j along the rotor drift direction (red solid line in voltage snapshot). Insert is a magnified view showing that h*J close to the rotor core is larger on the LA side than on the PV side (red arrows). (B) Left: snapshots of voltage (gray scale) and h·j (color coded) maps at two instances: when the rotor wave front near the core is directed from PV to LA (top) and about half cycle later and when it is directed from LA to PV (bottom). Blue dots represent the tip of the rotor and white dots are 4 mm from tip, toward LA (top) and toward PV (bottom). Right: Single pixel time series showing simultaneous voltage and h·j at the white dots. Red arrows indicate the times of the voltage and h·j snapshots on the left side of the panel (reproduced from Ref. 19).

Figure 5B shows snapshots of the voltage and h·j at a moment when the wave front near the rotor tip is propagating toward the LA (top) and a half cycle later, when that wave front is propagating toward the PV (bottom). The h·j snapshots clearly demonstrate that the wave propagating toward the LA is facing a higher h·j (red) as compared with the wave propagating toward the PV. The h·j gradient during the drift is further confirmed by plotting the time course of the voltage and h·j at 2 pixels flanking momentarily the SP of the rotor: one on the LA side and the other on the PV side, 4 mm away from the tip. Those dynamical plots show that the peak h·j for each cycle at the LA side is always >0.4 and at the PV side is always <0.4. As illustrated, the alternating SP drift is larger during propagation toward the PV (low h·j) than during its propagation toward the LA (high h·j). The result is a net drift toward the lower h·j, ie, less excitable, region. Importantly, it is noticeable that the presence of heterogeneity at the core area is critical for the rotor drift; in the absence of heterogeneity at the core, a rotor may be stable. ¹⁹

Figure 6 summarizes the relationship between rotor drift and spatial distribution of excitability factors during reentry in the three conditions simulated. Figure 6A shows for illustration purposes a map with the distribution of the time-averaged h·j peak values for the model with greatest excitability gradient (Condition III) and shows that those values are about twofold larger in the LA edge as compared with those at the PV edge (0.8 and 0.4, respectively). The graph in Figure 6A presents the spatial profiles of the h·j and the drift direction for the three conditions. It is seen that the two conditions with lowest density of I_K1 at the PV edge (blue and red) have profiles with reduced h·j at that edge as well, in contrast to the condition without the I_K1 gradient, where h·j is maximal at the PV edge. Overall, the directionality of the drift in the three conditions as indicated by the superimposed arrows is fully consistent with the rotor attraction by any region with lowest h·j values where the steepest h·j gradient (red) corresponds to the fastest drifting rotor, as indicated in Figure 4B (red).

Figure 6.

Gradients of sodium availability and AP measurements during rotor drift predict drift direction. (A) Left: map of averaged (h·j)_peak (<(h·j)_peak>) during rotor activity under Condition III (only I_K1) on a flat model is showing a LA to PV decrease. Right: PV-LAJ spatial profiles of <(h·j)_peak> for the three conditions modeled. Superimposed arrows indicate direction of rotor drift for each condition, consistently toward the region of lowest <(h·j)_peak> (lowest sodium availability). (B) AP measurements of excitability during reentry. Left: A voltage snapshot and APs from the LA (red) and PV (blue) edges (asterisks on map) for the three conditions simulated. Right: upstroke velocity (dV/dt_max) and MDP at the LA and PV region with direction of drift (arrows) for the three conditions. Drift is consistently toward lowest upstroke velocity and most positive MDP regions (reproduced from Ref. 19).

Next, Figure 6B quantifies metrics associated with the AP measurements of excitability: we focus on the dV/dt _max and MDP, which plays a role in determining the availability of I_Na during the membrane depolarization. On the left side, we present samples of aligned pairs of APs recorded in locations (asterisks) near the LA (red) and PV (blue) edges during reentry in the three conditions. As can be appreciated, each condition presents a distinct heterogeneity in its APs as quantified on the right side of the panel, but for the three conditions, the drift direction is toward the regions with lowest excitability as determined by the slowest dV/dt_max and most positive MDP. ²²

Figures 3 and 4 demonstrated the important role of I_K1 in determining the direction of the rotor drift in the PV-LAJ. In Figure 7 we study the effect of various relevant current–voltage relationships of I_K1 on such drift direction. Figure 7A shows four relationships between the current density and transmembrane voltage (I–V relationship) for different I_K1 levels, including up- and downregulation, as well as their corresponding APs showing different APDs and MDPs. The four I–V relationships were incorporated in a PV-LAJ model with Condition I, and rotor activity as well as AP parameters were tracked. As can be observed from Figure 7B, increasing or decreasing the I_K1 increases or decreases the average rotor frequency, respectively. Figure 7C demonstrates that h·j profiles across the PV-LAJ vary in both levels and gradient directions as a consequence of altering the I–V relationships. In particular, it is noticeable that the rotor drifts toward the PV for a broad range of I_K1 levels, whereas a reversal in the drift and h·j gradient direction occurs only at a 75% reduction in I_K1 (a behavior similar to that of a uniform I_K1 shown under Condition II in Fig. 4).

Figure 7.

Sensitivity of the rotor drift direction to I_K1 I–V profile. (A) I_K1 I–V curves for normal (control) and other levels of I_K1 density. Steady-state APs induced at 2 Hz are illustrated. (B) Rotation frequency in LA homogenous model for each I_K1 density level. (C) Average peak I_Na availability profiles show drift toward lower excitability (arrows) in a flat model with Condition I. (D) Analysis of MDP and dV/dt_max for each profile at single pixels in the LA and PV (4 mm from boundaries) under Condition I (reproduced from Ref. 19).

Up- and downregulation of I_K1 also hyperpolarize and depolarize MDP as well as increase and decrease upstroke speed, respectively (Fig. 4), with consistent drifts toward higher MDP and lower upstroke velocity. Overall, the consequences of altering the I–V properties of I_K1 in the heterogeneous models on MDP and upstroke velocity are complex and suggest a nonmonotonic dependency. Nevertheless, further alterations in extracellular potassium concentration in the range of 4-7 mM, as well as alteration in the I–V profile to resemble closely those of Kir2.1 or Kir2.3, affect excitability and rotation frequency but do not alter direction of h·j gradients and drift directions toward the PV. ¹⁹

Discussion

The findings presented here suggest the following: (1) the heterogeneous distribution of transmembrane currents in the PV-LAJ plays a major role in the preferential localization of rotors near or at the PVs, and the spatial dispersion of I_K1 is particularly important in determining such attraction. (2) Rotor drift in the PV-LAJ can be attributed to an excitability (maximum I_Na availability) gradient in the vicinity of the rotor pivoting point.

Mechanisms of AF and Rotor Dynamics

Our understanding of AF in individual patients would benefit from knowledge of how driving rotors form and then become stable or unstable, under the conditions of multifactorial substrate heterogeneity.^25,26 Paroxysmal AF in patients and in isolated normal sheep hearts has been found often to depend on fast rotors localized to the posterior wall of the LA and the PVs region with fibrillatory conduction toward the rest of the atria.^1,2 However, the ionic properties underlying those rotors' formation and drift remain unclear.^7,27 Moreover, recent studies have indicated that in some patients, rotors that drive the AF may reside outside of the PV area. ⁴ Our simulations here demonstrate that the spatial distribution of ionic currents found in the canine PV-LAJ is conducive to attracting rotors to the PV region. In addition, we also demonstrate that this attraction toward the LA can be reversed or arrested, if certain ionic currents are altered, which in turn may explain the variability in the location of rotors found in different patients. However, AF may involve several coexisting rotors at any given moment. In these cases, in addition to the drift imposed by the underlying substrate, the faster rotors can also exert an overriding influence on the slower rotors, ²⁸ and the combination of these two factors on rotor dynamics warrants further investigation.

Substrate Heterogeneity and Rotor Drift

In addition to the role of restitution characteristics of cells in rotor stability, ²⁹ simulations with gradient of excitability showed spiral wave drifting in the direction of the region exhibiting lower excitability and velocity, with additional perpendicular component depending on the rotor chirality, ⁸ excitability, and repolarization, ⁹ regardless of the details of the initial conditions. ³⁰ Simulations using more biophysically detailed ionic models found that for a fixed gradient in APD imposed by linearly varying potassium currents, the velocity of the drift of a rotor is a function of the magnitude of the gradient ¹⁰ and a steep gradient in APD can lead to conduction block of premature beats. ³¹ Our simulations in a flat and cylindrical (to account for a drift not fully aligned with the gradient direction) models ¹⁹ are in agreement with all those previous studies but refine the previous prediction that rotors would drift toward regions with longer APD^9,32,33 to only >80% repolarization measured at frequencies close to the rotor frequency. Our simulations also show (Figs. 6 and 7) the drift predictability of various measures of excitability, which includes MDP, dV/dt_max, and h·j (I_Na availability).

I_K1, I_Kr, and Rotor Dynamics During Fibrillation

I_K1 and I_Kr density gradients in the dog PV-LAJ are found to be opposing each other. ¹⁴ As recent studies have shown that these two currents are important in rotor dynamics and AF,^{12,13,34–37} the ionic mechanisms leading to the propensity of the PV region to favor rotor activity^16–18 became complex. We show for the first time ¹⁹ a clear propensity of the currents distribution in the PV-LAJ to attract rotors to the PVs and the dominant role for the I_K1 dispersion over all currents, and in particular I_Kr, in determining the localization of a rotor in that area and open the possibility that interplay between I_K1 and I_Kr may be important for the differential localization of rotors in AF.

Comparing with other recent studies, the drift toward low I_K1 in our study is fully consistent with the simulations by Kneller et al. ²¹ , Comtois et al. ³⁸ , and Comtois and Nattel ³⁹ who studied the effect of artificial heterogeneity in the inward rectifier I_K,ACh on AF dynamics. Their simulations also suggest that while rotors accelerate their rotation frequency with increasing I_KACh, ²⁷ the low I_K,ACh regions are the ones that attracts rotors. ²¹ Other studies also described rotor attractions to long APD regions, albeit with a sharp heterogeneities. ⁴⁰ A recent study by Sekar et al. ⁴¹ utilized circular monolayers with overexpression of I_K1 either in a central circular island or in its periphery to show that rotors stably pivot around the island regardless of the relative level of I_K1. In that study, the gradient, however, was very sharp relative to the size of the rotor core, and the preparation was highly symmetrical, which may explain why that study did not show a preferential anchoring of rotors to either low or high inward rectifying K⁺ current levels as observed in this and other studies.^21,38,42 Finally, in a recent study in cardiomyocytes monolayers with heterogeneous I_Kr expression, stable rotors localized to the region with highest expression of I_Kr. ³⁷ Those stable rotors did not drift as in our simulations since they reside in a relatively uniform large region, in accordance with our simulations presented elsewhere. ¹⁹

Limitations

We study a specific set of membrane kinetic models (CRN-K) with a Boltzmann distribution of the current densities across the PV-LAJ preferred for extensive validation of propagation properties^22,43 over a more recent and detailed model requiring adjustments.^44,45 Experimental or clinical data on ionic properties and dispersion in the atria are scarce; we focused on the effect of reported ionic currents data for the dog. However, attraction or repulsion of rotors by the PVs at the PV-LAJ may be affected by factors other than those studied here. For example, the heterogeneity in the intrinsic cellular properties may have different effective heterogeneity in refractory and excitability ³⁰ depending on structural intercellular coupling, ⁴⁶ fibrosis, ⁴⁷ or the size of the medium.^31,39,48,49 Further, the drift of rotors may be influenced by accumulations of intra- or extracellular ions, as has been shown to occur in AF. ⁵⁰ Our study ruled out that PV narrow funnel-like anatomy reverses the ionic-induced attraction to the PVs ¹⁹ ; however, additional anatomical factors such as wall thickness,^51,52 the fiber bundles, ⁵³ or fibrosis ⁷ may also regulate the drift of rotors, possibly even counteracting the drift trend caused by the ionic gradients. To mitigate these limitations, we focused in our study on conditions relevant only to paroxysmal AF, prior to any remodeling and fibrosis, and incorporated various possible I–V relationships to substantiate our conclusions regarding the I_K1 dominance and drift prediction. ¹⁹ Our study nevertheless should be considered only as a first step in elucidating the concept of heterogeneity-induced drift and needs to be tested in future experimental studies.

Conclusions

Consistent with experimental and clinical studies on paroxysmal AF, simulations in an ionically heterogeneous model of the PV-LAJ showed rotor attraction toward the PVs. Our simulations suggest that I_K1 heterogeneity across the PV-LAJ is dominant compared to other currents in conveying the drift direction through its effect on refractoriness and excitability.

Author Contributions

Conceived and designed the experiments: OB. Wrote the first draft of the manuscript: OB. Developed the structure and arguments for the paper: OB. Made critical revisions and approved final version: OB. The author reviewed and approved of the final manuscript.

Abbreviations

AP

action potential

APD

action potential duration

AF

atrial fibrillation

Condition I

all currents varied spatially in the model

Condition II

all currents, except I_K1, varied spatially in the model

Condition III

only I_K1 varied spatially in the model

CRN-K

Courtemanche-Ramirez-Nattel and Kneller

CV

conduction velocity

DF

dominant frequency

dV/dt_max

maximum upstroke velocity

G _x

maximum conductance of channel x

h·j

the fraction of sodium channels available for activation

h·j _peak

the maximum fraction of sodium channels available for activation during a cycle

I _K1

inward rectifier potassium current

I _K,ACh

acetylcholine modulated inward rectifier potassium current

I _Kr

rapid delayed rectifier potassium current

I _Ks

slow delayed rectifier potassium current

I _to

outward transient potassium current

I _Na

sodium current

I _CaL

L-type calcium current

LA

left atrium

LV

left ventricle

MDP

minimum diastolic potential

PVs

pulmonary veins

PV-LAJ

LA and PV junction

RV

right ventricle

SP

singularity point. (Used to indicate the pivoting location of a rotor)

S1, S2

two sequential stimulations at particular times and locations

TSP

time–space plot.