Individual Contributions of Cardiac Ion Channels on Atrial Repolarization and Reentrant Waves: A Multiscale In-Silico Study

The excitation, contraction, and relaxation of an atrial cardiomyocyte are maintained by the activation and inactivation of numerous cardiac ion channels. Their collaborative efforts cause time-dependent changes of membrane potential, generating an action potential (AP), which is a surrogate marker of atrial arrhythmias. Recently, computational models of atrial electrophysiology emerged as a modality to investigate arrhythmia mechanisms and to predict the outcome of antiarrhythmic therapies. However, the individual contribution of atrial ion channels on atrial action potential and reentrant arrhythmia is not yet fully understood. Thus, in this multiscale in-silico study, perturbations of individual atrial ionic currents (INa, Ito, ICaL, IKur, IKr, IKs, IK1, INCX and INaK) in two in-silico models of human atrial cardiomyocyte (i.e., Courtemanche-1998 and Grandi-2011) were performed at both cellular and tissue levels. The results show that the inhibition of ICaL and INCX resulted in AP shortening, while the inhibition of IKur, IKr, IKs, IK1 and INaK prolonged AP duration (APD). Particularly, in-silico perturbations (inhibition and upregulation) of IKr and IKs only minorly affected atrial repolarization in the Grandi model. In contrast, in the Courtemanche model, the inhibition of IKr and IKs significantly prolonged APD and vice versa. Additionally, a 50% reduction of Ito density abbreviated APD in the Courtemanche model, while the same perturbation prolonged APD in the Grandi model. Similarly, a strong model dependence was also observed at tissue scale, with an observable IK1-mediated reentry stabilizing effect in the Courtemanche model but not in the Grandi atrial model. Moreover, the Grandi model was highly sensitive to a change on intracellular Ca2+ concentration, promoting a repolarization failure in ICaL upregulation above 150% and facilitating reentrant spiral waves stabilization by ICaL inhibition. Finally, by incorporating the previously published atrial fibrillation (AF)-associated ionic remodeling in the Courtemanche atrial model, in-silico modeling revealed the antiarrhythmic effect of IKr inhibition in both acute and chronic settings. Overall, our multiscale computational study highlights the strong model-dependent effects of ionic perturbations which could affect the model’s accuracy, interpretability, and prediction. This observation also suggests the need for a careful selection of in-silico models of atrial electrophysiology to achieve specific research aims.


Introduction
The cardiac action potential (AP) results from a complex dynamic behavior of ionic currents within a cardiomyocyte. Together, they maintain the excitation, contraction, and relaxation of a cardiomyocyte. In a physiological condition, the activation of voltage-gated Na + channels induces the opening of other voltage-gated ion channels, including the L-type Ca 2+ channels, allowing access of Ca 2+ from the extracellular space to the cytoplasm. Subsequently, the Ca 2+ 2 of 11 influx initiates Ca 2+ -induced Ca 2+ release (CICR), releasing more Ca 2+ from the sarcoplasmic reticulum to the cytosolic space to exert their function in cardiomyocyte contraction and other Ca 2+ -dependent signaling processes [1]. Such an activation of fast Na + channels is denoted as a positive deflection (rapid depolarization; phase 0) of the cardiac AP, whereas the opening of the L-type Ca 2+ channels modulate the plateau phase (phase 2) of the AP. During repolarization, several K + channels, including the transient-outward and delayed-and inward-rectifier K + channels, are activated. In the atria, some atrial-specific K + channels (e.g., Ca 2+ -activated, two-pore domain and ultra-rapid delayed-rectifier K + channels) have also been reported to hold an important role in the AP repolarization. Such a repolarization is displayed as a negative deflection of the AP (phase 1 and 3), restoring the membrane potential to the resting state (phase 4). Meanwhile, in the presence of disease-associated ionic remodeling, the collaborative efforts between ion channels could be disrupted, allowing cardiac arrhythmias to occur and persist. Therefore, AP morphology and AP duration (APD) are commonly considered as cellular markers of cardiac arrhythmias.
Computational modeling of atrial electrophysiology has been employed to study the mechanisms of atrial arrhythmias, to predict the outcome of antiarrhythmic interventions, and to assist decision-making process in the clinic [2,3]. For example, multiscale computational modeling was employed to study the ethanol-associated atrial arrhythmogenesis [4] and to support individualized planning for catheter ablation procedures of atrial arrhythmias [5]. However, the distinct formulations within each atrial model could affect the simulated results and therefore might influence the accuracy of in-silico models to predict the outcome of specific interventions. Moreover, a better understanding on the individual contribution of atrial ion channels on atrial repolarization is a prerequisite for designing an appropriate channel-targeted antiarrhythmic therapy.
Therefore, this multiscale in-silico study sought to explore the individual contribution of atrial ionic currents on AP repolarization and reentrant waves behavior in two distinct insilico models of atrial cellular electrophysiology, and to test several potential antiarrhythmic strategies to destabilize AF-associated reentrant spiral waves in silico.

Zero-Dimensional Computational Modeling
The individual contribution of nine major ionic currents (fast Na + current [I Na ], transient-outward K + current [I to ], L-type Ca 2+ current [I CaL ], ultra-rapid delayed-rectifier K + current [I Kur ], rapid delayed-rectifier K + current [I Kr ], slow delayed-rectifier K + current [I Ks ], inward-rectifier K + current [I K1 ], Na + -Ca 2+ exchange current [I NCX ] and Na + -K + pump current [I NaK ]) were assessed in two in-silico human atrial cardiomyocyte models, namely Courtemanche et al. [6] and Grandi et al. [7]. To simulate the cellular effect of AF-associated ionic remodeling, the AF variant of the Courtemanche model, which was made by incorporating the chronic-AF-associated electrical remodeling of transmembrane currents reported by Grandi et al. [7] (i.e., I Na −10%, I to −80%, I CaL −50%, I Kur −55%, I Ks +100%, I K1 +100% and I NCX +40%) were employed. AP simulations were performed at a 1 Hz pacing frequency (basic cycle length [BCL] of 1000 ms) and quasi-steady-state APD for the baseline models was obtained following 100 beats of pre-pacing. All original models were obtained from CellML and all simulations were performed in Myokit [8].

Two-Dimensional Computational Modeling
Reentrant spiral waves were simulated using an S 1 S 2 induction protocol in homogeneous tissue of 4 × 4 cm (200 × 200 units) following 100 beats of 1 Hz pre-pacing at the cellular level. The tissue CV was maintained around 45-50 cm/s for the baseline models and allowed to vary in the presence of underlying ionic perturbations. The first stimulus (S 1 ) was initiated from left to right to generate a normal excitation wave. Subsequently, the second stimulus (S 2 ) was applied to the upper-left quadrant of the tissue, generating an additional wavefront that can interact with the tail of the preceding wave, producing reentry in a vulnerable substrate [4]. The vulnerable window was evaluated to assess both the inducibility and stability of reentrant arrhythmias under different conditions. The size of the vulnerable window indicates the inducibility of reentry, whereas the duration of reentry within the vulnerable window indicates the stability of reentrant arrhythmias. A stable reentry was defined as a reentrant spiral wave that persisted after 12 s.

Individual Contribution of Ion Channels on Atrial Action Potential Repolarization
To explore the contribution of individual atrial ionic currents on AP repolarization, a sensitivity analysis was conducted by varying the maximum conductance (G max ) of ionic currents of interest up to 400% of the default values (Figures 1 and 2). Figure 1 display sensitivity plots describing the changes of APD at 90% repolarization (APD 90 ) due to ionic current alterations (0% to 400%). In the Courtemanche model ( Figures 1A and 2A), all currents, except I Na , had a meaningful effect on the AP repolarization phase, whereas in the Grandi human atrial cardiomyocyte model ( Figures 1B and 2B), I Kr and I Ks had limited contribution on AP repolarization. Notably, a 50% reduction of I to abbreviated the APD 90 in the Courtemanche model (268 ms vs. 293 ms [−8.5%]), whereas APD 90 was prolonged (331 ms vs. 302 ms [+9.6%]) in the Grandi model. The Grandi model also exhibited a high sensitivity to I CaL perturbations, with I CaL above 150% of the baseline model resulting in plateau arrest ( Figure 2B). Meanwhile, the alteration of I Kur affected the entire phase of repolarization in the Grandi model, but its effect was primarily in the early phase of repolarization in the Courtemanche model. In both atrial models, I K1 perturbations notably affected the late phase of repolarization and resting membrane potential (RMP), strongly modulating APD 90 . In the Courtemanche model, perturbations in I NCX and I NaK that prolonged APD during the early repolarization phase shortened APD 90 .   [6] and Grandi et al. [7]] were used to assess the contribution of nine major ionic currents (i.e., I Na , I to , I CaL , I Kur , I Kr , I Ks , I K1 , I NCX and I NaK ) on atrial repolarization. Channel maximum conductance (G max ) were scaled from 50% to 400% of the default value.

The Effects of Ionic Current Conductance Perturbation on Reentrant Spiral Waves
The effect of an antiarrhythmic therapy in the cellular level could be different than that in the tissue and organ levels due to the potential influence of cell-cell interactions. Therefore, to improve our understanding on the consequences of ion-channel perturbations on the behavior of reentrant arrhythmias, we performed sensitivity analyses in homogenous atrial tissues, assessing the individual contribution of nine major ionic currents.
Inhibitions of I Na shifted the vulnerable windows to larger S 1 S 2 intervals and widened the size of the vulnerable windows, indicating a higher inducibility of re-entry in both the Courtemanche and Grandi models ( Figure 3A,K). Similar to the cellular sensitivity analysis (Figure 2), the response to I to block was different in the Courtemanche and Grandi models.
In the Courtemanche tissue model ( Figure 3B), I to block produced a shift of the vulnerable windows to the earlier S 1 S 2 intervals, without affecting the reentry inducibility. In contrast, in the Grandi model ( Figure 3L), I to block shifted the vulnerable windows to larger S 1 S 2 intervals, consistent with the cellular APD response. The upregulation of I to up to 200% in the Courtemanche model reduced the inducibility of reentrant waves, whereas increasing I to to 300% markedly shifted the windows to earlier S 1 S 2 intervals. However, the stability of reentry was not affected by I to perturbations. Meanwhile, increasing I to in the Grandi model shifted the vulnerable windows to earlier S 1 S 2 intervals. At 200% I to , the stability of reentry increased, and a stable reentry was attained with 300% I to . Reductions in I CaL shifted the vulnerable windows to earlier S 1 S 2 intervals in both models, with some instances of stable reentry in the Grandi model following 25-50% I CaL (Figure 3M), whereas I CaL upregulation shifted the windows to larger S 1 S 2 intervals in the Courtemanche model ( Figure 3C) and to shorter intervals in the Grandi model ( Figure 3M). Of note, the I CaL upregulation cannot exceed 150% in the Grandi model due to repolarization failure at higher current densities. Alterations of I Kr and I Ks did not have a major effect on AP properties in the Grandi model (Figure 2). Similarly, at the tissue level, only I Kr upregulation shifted the vulnerable windows to earlier S 1 S 2 intervals, with no effect on the inducibility and stability of reentry ( Figure 3O). Meanwhile, in the Courtemanche model, I Kr upregulation (300%) stabilized reentrant waves ( Figure 3E). Similarly, the upregulation of I K1 could stabilize reentrant arrhythmias in the Courtemanche model ( Figure 3G), while no stable reentry was documented in the Grandi model ( Figure 3Q). Figure 3J shows the consequence of AF-related electrical remodeling incorporated in the Courtemanche atrial model on the behavior of atrial reentry. Figure 4 shows some examples of the effect of ion-channel perturbations on reentrant spiral waves. As indicated above, 300% I Kr , 300% I K1 and AF electrical remodeling triggered stable reentry in the Courtemanche model, whereas 300% I to and 25% I CaL resulted in similar behavior in the Grandi model. The rotor cores ( Figure 4A-H, right panels) meandered less in the stable reentry as compared to the unstable ones ( Figure 4A,G), consistent with previous work [9]. Overall, our tissue-scale sensitivity analyses, as previously presented [10], indicate a strong model dependence of the effect of ionic perturbations on reentrant arrhythmias, making it difficult to interpret their exact consequences in silico. These findings emphasize the importance of developing a personalized in-silico model to assess the patient-specific arrhythmogenic risk and predict the individual response to antiarrhythmic therapies.

In-Silico Assessment of Antiarrhythmic Properties of Ionic Current Perturbations
Taking into account the insight from the tissue sensitivity analyses, we tested the efficacy of three different modalities to treat AF-associated electrical remodeling in homogenous atrial tissues using the Courtemanche human atrial cardiomyocyte model: the recovery of AF-related I to remodeling, recovery of AF-related I CaL remodeling and 80% I Kr block, mimicking the effect of class III antiarrhythmic drugs (e.g., dofetilide and ibutilide). Reversing I to remodeling shifted the vulnerable windows to the earlier S 1 S 2 intervals and slightly reduced the inducibility of reentrant waves. However, some episodes of stable reentry were still detected ( Figure 5A). On the contrary, reversing I CaL remodeling shifted the vulnerable windows to the larger S 1 S 2 intervals and destabilized reentry, with no stable reentry remained. Similarly, 80% block of I Kr also countered the AF electrical remodeling and yielded similar effects to I CaL remodeling recovery. Based on these findings, we further explored whether I CaL recovery and I Kr block had identical antiarrhythmic properties in the virtual atrial tissue. To achieve this aim, we tested the efficacy of both interventions in two different circumstances, namely acute and chronic settings. In an acute setting, the pharmacological intervention is applied after the reentry is initiated. Meanwhile, in the chronic setting, the treatment modality is applied before the reentry occurs. Clinically, the former resembles an acute treatment of cardiac arrhythmias, whereas the latter mimics the ability of a treatment modality to prevent the reoccurrence or exacerbation of cardiac arrhythmias. Figure 5B depicted the representative snapshots of reentrant waves (S 1 S 2 interval = 230 ms) demonstrating that both modalities were effective in destabilizing reentrant arrhythmias in the chronic setting. However, Figure 5C,D clearly showed that only 80% I Kr block was effective to destabilize and terminate reentrant arrhythmias in an acute setting, whereas acute recovery of I CaL remodeling preserved the existence of the reentrant spiral wave.

Discussion
The repolarization of atrial action potential is maintained by numerous ionic currents [11][12][13][14][15][16][17]. In canine atrium preparations, lidocaine (a class IB antiarrhythmic drug that predominantly blocks I Na ) minimally shortened APD 90 but caused a significant effective refractory period (ERP) prolongation. Meanwhile, E-4031 (a specific I Kr blocker) prolonged both APD 90 and ERP, whereas the combination of the two displayed a synergistic AFsuppressing effect [18]. The effects of I Na and I Kr inhibition on atrial cardiomyocytes were also seen in our study, in which I Na inhibitions had negligible impact on cellular APD 90 ( Figure 2) but significantly shifted the vulnerable windows to the larger S 1 S 2 intervals (Figure 3), suggesting a prolongation of tissue ERP. Additionally, the reduction of I Kr extended the APD 90 and shifted the vulnerable windows to the larger S 1 S 2 intervals in the Courtemanche human atrial model. Several studies have also highlighted the importance of I to in the repolarization of atrial cardiomyocytes. A human right atrial experiment revealed that I to produced pronounced changes on early repolarization of atrial AP and force generation [19], and its contribution was effective at all physiological heart rates in humans [20]. Ni et al. [21] also explored the contributions of I to in silico and demonstrated that I to modulated the atrial APD restitution and its alteration contributed to APD alternans. The in-silico effect of I to upregulation on APD seemed to be consistent across computational models of atrial cardiomyocytes [17], including our current work, in which we observed reductions of APD in both early and late repolarization. At tissue scale, the role of I K1 in reentrant wave stabilization has been documented in both computational and experimental studies [22]. Computer simulations indicated that in addition to the APD shortening, I K1 upregulation also increased the availability of I Na , which further accelerated the rotor. Experimentally, a similar observation was documented in I K1 -overexpressed mice, in which the rotor was persistent and faster than the wildtype animals [22]. Consistent with our simulation in Figure 5D, in addition to its effect on the prevention of AF recurrence, dofetilide has been widely used as an effective rhythm-control strategy for acute AF cardioversion. Spontaneous conversion of persistent AF was commonly seen within three days of dofetilide administration [23]. Meanwhile, although currently no available pharmacological agent could revert the pre-existed AF-associated I CaL remodeling, in the future, a drug that upregulates the I CaL density, possibly through an alteration of channel's kinetics or expression, could be beneficial to prevent the AF recurrence, as we demonstrated previously in Figure 5A,B. Alternatively, gene therapy might also be useful to revert AF remodeling through several ways, as previously discussed by Liu and Donahue [24].
Next, our multiscale computational study also indicates strong model-dependent effects of ionic perturbations on atrial electrophysiology. For example, in the Grandi atrial model, the contributions of I Kr and I Ks on atrial repolarization and reentrant waves behavior were minor, affecting the suitability of this model for in-silico drug cardiotoxicity screening and the mechanistic investigation on the arrhythmic consequences of disease-associated I Kr and I Ks remodeling. Of note, experimentally, the presence of both rapid and slow components of delayed-rectifier K + (Kr and Ks) channels has been previously reported in human atrial cardiomyocytes isolated from right atrial appendages extracted at the time of coronary artery bypass surgery [25]. Additionally, in Figures 2 and 3, we also exemplified the model-dependent effect of I to inhibition on atrial AP and reentrant waves. In the Grandi human atrial cardiomyocyte model, 50% I to inhibition prolonged APD and shifted the vulnerable windows to larger S 1 S 2 intervals, whereas in the Courtemanche human atrial model, I to inhibition resulted in the shortening of APD 90 and the shifting of vulnerable windows to earlier S 1 S 2 intervals. Experimentally, the inhibition of I to by 4-aminopyridine (4-AP) in human atrial cardiomyocytes obtained from human right atrial appendages extracted at the time of bypass surgery also abbreviated the phase 1 of AP (increasing the plateau height) and shortened APD 90 without any effect on RMP [19], consistent with our observation in the Courtemanche human atrial model ( Figure 1A).
Overall, the agreement of our findings with previous experimental studies demonstrates that in-silico modeling is a powerful and robust tool to study the mechanistic background of atrial arrhythmias (e.g., to explore the effects of ion-channel perturbation on atrial AP and reentrant arrhythmias) and to guide arrhythmia therapy in the clinic (e.g., to assess the acute and chronic effects of anti-AF medications). Nonetheless, a careful assessment on the appropriateness of each mathematical model to do specific tasks is warranted due to the strong model dependence. Norbert Wiener's infamous quote in 1945-"The best material model for a cat is another, or preferably the same cat" highlights such a basic limitation of computational models that no in-silico model could perfectly match the modeled biological systems (e.g., cells, tissues or organ). Also, it is important to note that an in-silico model is commonly validated to a certain extent and is designed to fulfil specific research objectives, therefore the application beyond those circumstances (i.e., the validated range of operation) requires careful consideration [26].
Finally, several limitations of this study need to be acknowledged. First, when simulating the AF-related remodeling, we did not incorporate the AF-associated remodeling of Ca 2+ -handling proteins, as well as AF-induced structural remodeling. Second, we also did not incorporate atrial-specific ion channels in the models since they are not available in the original models and for acetylcholine-activated inward-rectifying K + currents, since we did not simulate the effects of parasympathetic/vagal response in this study. Third, we performed the tissue simulations in a small sized tissue (4 × 4 cm) due to limited computational power. Our preliminary assessment (Supplemental Figure S1) revealed that tissue size could influence the stability of reentry, presumably due to the availability of excitable area for wave propagation, allowing the meandering rotor to persist. However, we speculate that this issue would not affect the overall results of this study since all of the comparisons were made in the same tissue size, although the confirmation of this notion in a larger tissue size is warranted. Additionally, in the future, an extension to organ-level simulations could be performed.

Conclusions
Atrial repolarization is regulated by multiple ion channels, which further modulate the behavior of reentrant spiral waves in the tissue level. Our multiscale computational study highlights the strong model-dependent effects of ionic perturbations which could affect the model's accuracy, interpretability, and prediction. This observation also suggests the need for a careful selection of in-silico models of atrial electrophysiology to achieve specific research aims.

Data Availability Statement:
Computer codes used in this study, such as the Myokit files for the zeroand two-dimensional simulations, are available from: https://github.com/henrysutanto (accessed on 2 December 2021).

Conflicts of Interest:
The authors declare no conflict of interest.