To obtain data for the population of models, samples from the right atrial appendages from 149 patients diagnosed with chronic AF in which antiarrhythmic medication was interrupted before the study were available (Wettwer et al., 2013; Sánchez et al., 2014).

Briefly, patch-clamp was performed in all the samples obtaining the values for different currents. A total of six biomarkers were used to quantify variability in action potentials (AP) including AP duration at 20, 50, and 90% of repolarization (APD20, APD50, and APD90, respectively), AP amplitude (APA), resting membrane potential (RMP), and AP plateau potential at 20% of APD20 (V20). The maximum and minimum values of these biomarkers at a pacing frequency of 1 Hz are presented in Supplementary Table 1.

To build the computational population of human AF models, different combinations of ionic currents were generated from the experimental data described earlier. In further detail, Latin Hypercubic Sampling (LHS) (McKay et al., 1979) and the baseline AF model developed by Koivumaki et al. were used for its generation.

A total of nine parameters were varied from –50 to +100% of their original value: fast sodium (Na⁺) ionic conductance (g_Na), sodium-potassium ion (Na⁺–K⁺) pump (I_NaK), inward rectified K⁺ current (gK1), L-type calcium ionic conductance (g_CaL), ultrarapid outward ionic conductance (g_Kur), Ca²⁺-dependent K⁺ current (I_KCa), and Na⁺ and K⁺ extracellular concentration and the diffusion coefficient of the reaction-diffusion equation (Supplementary Figure 2). LHS produced 500 different combinations of these nine parameters from the initial 149 patients. Simulations for these 500 profiles were calculated on planes of 8 × 256 nodes to evaluate the AP metrics. The models were simulated by pacing at 1 Hz (using a 3 ms stimulus duration, twice diastolic threshold amplitude). The APs of three cells along the plane (cells 500, 620, and 748) were analyzed following a train of 15 periodic stimuli, and the last 5 periodic stimuli were considered in order to ensure a steady-state. Only the models that fitted into the experimental constraints (Supplementary Table 1) were considered representative human electrophysiological models for the final population, resulting in 127 final profiles (Simon et al., 2017).

Different atrial models have been described at a unicellular level to characterize the electrophysiological response. In this study, we implemented the Koivumaki model with Skibsbye modifications (Skibsbye et al., 2016) that included a reformulation in sodium current to characterize sodium channel inactivation, adjusted transient outward potassium current and L-type Ca⁺ current, and included the small conductance calcium-activated potassium current (I_KCA). In addition, the model used for simulations included AF remodeling, achieved by modifying the following ionic currents: L-type Ca²⁺ (I_CaL) decreased by 55%, transient outward current (Ito) decreased by 62%, rapid delayed rectifier potassium channel (I_Kur) decreased by 38%, inward-rectifier potassium channel (I_K1) increased by 62%, Na/Ca exchange current (I_NCX) increased by 50%, expression of sarcoplasmic/endoplasmic reticulum Ca²⁺-ATPase (SERCA) decreased by 16%, phospholamban to SERCA increased by 18% and sarcopilin to SERCA decreased by 40% as described in Koivumaki et al. (2014). This model was used for both the calibration of the population of models previously described and the rest of the experiments in the study.

Simulations were performed on 2D planes mimicking a sheet of cardiac tissue. Two different tissue sizes were implemented for this study: one square plane corresponding to a normal atrium (16 cm², 400 × 400 nodes) and another square plane corresponding to a dilated atrium (20.25 cm², 450 × 450 nodes) (Kou et al., 2014).

To connect the cells within the plane, the monodomain reaction-diffusion equation was implemented, assuming that tissue behaves as a functional syncytium where membrane voltage is propagated smoothly (Clayton and Panfilov, 2008):

where the ∇ corresponds to the gradient operator and D a diffusion coefficient with units distance² time^–1. By using this monodomain simplification, the tissue is considered to have an unlimited extracellular medium, so the extracellular resistivity can be neglected. The extracellular medium is isopotential and equal to zero for simplicity. Consequently, the membrane potential is the same as the intracellular potential. Planes were fully connected as shown in Figure 1, not including structures such as the pulmonary veins. The value of the diffusion constant, referred to as D in the above equation, was varied among the different profiles as part of the variability included in the population of models.

Both antiarrhythmic and proarrhythmic drugs were evaluated in the electrophysiological population of models in order to characterize the effect according to tissue size and ionic currents. Three different drugs were studied presenting different mechanisms and effects, namely, verapamil, flecainide, and isoproterenol. Briefly, verapamil is an antiarrhythmic drug that acts as a calcium blocker, flecainide is an antiarrhythmic drug that acts as a sodium blocker and isoproterenol is a proarrhythmic agent (β receptor agonist) that increases intracellular calcium. All three drugs are currently used in clinical practice and their effect, although established as proarrhythmic or antiarrhythmic, can vary from patient to patient (Bassett et al., 1997). A simple pore channel equation was used for drug implementation at the computational level including three channels for each drug that was modeled according to experimental IC50 values to calculate the block of the channel (Dempsey et al., 2014). The model used is described as follows:

where the G₀ represents the initial conductance of the channel for each of the profiles in the population of modes, [D_i] corresponds to the concentration of the drug and IC₅₀ is the concentration of the drug that reduces by 50% the channel current. For verapamil and flecainide, values from Crumb et al. (2016) were implemented. Isoproterenol was modeled to increase the permeability of the calcium current as stated in Vescovo et al. (1989). The numerical values for each modeled drug can be observed in Table 2, including the concentration at which the drug was modeled.

Simulations were performed implementing differential equations computed with a time step of 1 microsecond for the Euler method using in-house software written in C++ with CUDA parallelization and solved with a NVIDIA TESLA C2057 GPU, NVIDIA Corporation, Santa Clara, CA, United States. Rush Larsen method was developed for cell models as it offers stability to the problem by calculating the exact solution for the gating variables. Since all the equations governing the gating variables have a similar structure, the method uses the following expression to solve them:

Consequently, all equations for gating variables were solved by the previous expression whereas the rest of the equations were solved by the forward Euler Method.

Planes were simulated for a total of three impulses (S1) at 1 Hz followed by a fourth one (S2) for arrhythmia initiation. S2 was induced in the cells of the inferior left section of the plane, producing reentry in part of the models. For a given combination of ionic channels in the population of models, arrhythmia was induced when >1 complete rotational activity followed the S1–S2 induction protocol. In addition, rotor tracking was studied for two specific profiles that will be discussed in the “Result” section. The methodology for rotor tracking has already been described in previous publications of the group (Rodrigo et al., 2014). Briefly, phase maps of the simulations were calculated by using the Hilbert transform from which singularity points were calculated. A singularity point (SP) was defined as the point in a phase map that is surrounded by phases from 0 to 2π. Only those singularity points that were present for the duration of at least 1 full rotation were considered, as described in Rodrigo et al. (2014). Rotor track was defined as the connection between SPs across spherical layers at a given time. Only filaments that completed at least 1 rotation on the outermost surface were considered.

A total of 1,016 simulations were computed in this study corresponding to all different combinations of the population of models (127 profiles) simulated in different tissue sizes (two tissue sizes) and four different conditions (baseline conditions, two antiarrhythmic drugs, and one proarrhythmic drug), creating a database with different ionic conductance combinations and tissue size.

Random Forest, which is a decision algorithm consisting of a multitude of decision trees at training time, was implemented to output two possible outcomes: induced or non-inducible AF. This algorithm was trained including the eight variables of the population of models as an input and the presence of AF form simulations as an output to evaluate patterns that may lead to AF maintenance. The algorithm was trained with 80% of the data (N = 832) and 20% of the data was used for testing with a 5k-fold cross-validation.

From the complete population consisting of 127 different electrophysiological profiles, 80 profiles maintained the reentrant activity at baseline conditions in the normal tissue size and 87 in the dilated atrium. Complete quantification of the profiles maintaining reentry can be observed in Table 1, including the effect of the three simulated drugs. As shown, dilated tissue increased the number of profiles maintaining reentry in all cases, independently of the presence and type of drug.

The analysis of the simulations resulted in the identification of profiles that responded differently to the arrhythmia induction under the effect of the drug, as can be observed in arrhythmia induction Figure 2. Non-inducibility of the arrhythmia was observed in the majority of the profiles when the antiarrhythmic compounds, verapamil, and flecainide were used. However, in some of the profiles that did not induce the arrhythmia at baseline, the addition of one of those drugs induced rotational activity (Figure 2B). Specifically, the addition of flecainide gave rise to non-inducibility of the arrhythmia in 51 profiles and induced AF in seven profiles that were non-inducible at baseline. For verapamil, 53 and 8 profiles were non-inducible and induced the arrhythmia, respectively, in normal tissue size samples. Interestingly, all the profiles that showed proarrhythmic and antiarrhythmic effects in the flecainide scenario presented the same behavior in the verapamil scenario. Although verapamil and flecainide showed similar results, the overall action potential morphology was significantly different for the same profile under the effect of these drugs, as it can be observed in Figures 3, 4, where the curvature of the action potential is modified due to the effect of the drug at the ionic level. Isoproterenol showed a proarrhythmic effect increasing the number of AF maintaining profiles in eight cases.

Antagonistic effects of both antiarrhythmic and proarrhythmic drugs were observed on the population of models. Figure 3 presents a specific ionic profile simulated for all four conditions (basal and three drugs) in which reentry was maintained over time in the basal scenario. When the simulation was repeated under the effect of verapamil or flecainide, the arrhythmia was not induced. With the addition of isoproterenol, the arrhythmia was, not only induced but rotational activity presented a higher activation frequency.

For another specific profile shown in Figure 4, antagonistic effects were observed as follows: for a profile in which arrhythmia was not induced under basal conditions, verapamil and flecainide showed a proarrhythmic effect, meaning that the arrhythmia was induced, while isoproterenol did not induce the reentry for the complete simulation. Furthermore, for these aforementioned profiles, rotor tracking shows the lower area and higher complexity of the simulations in which the arrhythmia was maintained over time, exhibiting the stability of the reentry. Flecainide and verapamil terminated AF for the majority of the profiles in which AF was maintained under basal conditions, whereas isoproterenol has induced AF in new profiles which were not inducible in basal conditions. Interestingly, a small proportion of models presented AF induction under the effect of verapamil and flecainide (Figure 2B) despite the absence of AF induction at baseline.

We generated a final database of 1,016 simulations that included the basal and drug administration state for both tissue sizes. Therefore, these data represented: (1) the variables of the population of models, (2) the electrophysiologic change conferred by the drug administration; (3) the changes between normal and dilated tissue size. Then, the database was processed in order to train and calibrate a decision algorithm for drug effect prediction. As a result, the Random Forest algorithm, shown in Figure 5, was obtained. The algorithm had in total seven different consecutive layers, shown in Figure 5A that cluster similar profiles together for prediction of AF induction based on the conductance values in the form of a sunburst diagram. In this diagram, each level is represented by a concentric circle containing one or more variables from the population of models, where each level presents threshold values for decision making. To evaluate a specific profile, the decision algorithm starts the clustering process from the most inner circle, that corresponds to the G_K1 variable. For each variable, a threshold value has been defined that has to be compared with the value of the variable for the specific profile that is being evaluated. Specific values for each threshold can be consulted in the Supplementary Figure 1. For example, in the case of G_K1, this value corresponded to 9.12%. If the value of the profile is higher than 9.12%, the next level to be considered will correspond to the right part of the diagram and the next variable (High G_K1) to continue the characterization, corresponding to the diffusion variable. In contrast, if the value of the profile is lower than 9.12% for G_K1, the next level to be considered will correspond to the left part of the diagram and the next variable, identified as K₀. This process should be completed until the last layer or circle of the diagram, finding the path or cluster to which a given profile belongs. Figure 5B shows the paths or clusters that have been identified as inducible AF by the algorithm.

Interestingly, if the paths of clusters from the diagram are analyzed, it can be observed that one of the main contributions that lead to AF inducibility is the combination of profiles with high conductance of the inward rectifier K⁺ channel (K1), low diffusion, high concentration of extracellular potassium, high conductance for the sodium channel and low conductance of the ultra-rapidly activating delayed rectifier current (I_Kur). Another combination that led to AF was low conductance value for I_K1, low concentration of extracellular potassium, high conductance for the sodium-potassium pump, high conductance for the slow calcium channel, high conductance for the calcium-potassium pump, and high conductance for the I_Na channel.

Figure 6 shows an example of how the change in conductance of specific channels can affect the final path or cluster to which a specific profile belongs, altering the final outcome of the decision tree algorithm. In this figure, specific pathways of the decision tree are exemplified in each of the panels, highlighting how the change due to the addition of the drug affects specific levels that, depending on the final permeability of the channel, may be fundamental for the induction and maintenance of the rotational activity. Panels A, C highlight changes in calcium channel conductance, that can be affected by the addition of verapamil or isoproterenol, as both drugs were modeled to respectively, decrease or increase the permeability of this channel. Specifically, Panel A exemplifies how the increase in the permeability of the calcium channel, by the addition of isoproterenol, has a proarrhythmic effect. Conversely, a profile that presents high conductance of the calcium channel can be reduced by adding verapamil and changing to a state of non-inducibility of AF. Panel C shows another example in which the decrease of calcium permeability results in AF inducibility and the increase of calcium permeability increases the probability of AF non-inducibility.

Examples for flecainide, which affects sodium channel permeability, are shown in Panels B, D. Specifically, Panel B shows an antiarrhythmic effect of the addition of flecainide, by blocking sodium channels. Panel D shows a level of the decision tree algorithm at which, by changing sodium channel permeability, the probability of AF non-inducibility increases. Thus, the probability of inducing AF will be higher for those profiles with lower sodium permeability, showing a proarrhythmic effect of the drug.

The population of models has been widely used in the cardiac electrophysiology field for safety pharmacology, providing new platforms for the assessment of proarrhythmic effects of drugs (Passini et al., 2017; Kügler, 2020). In this study, we have analyzed how variations of the electrophysiological and anatomical characteristics can affect AF inducibility at the 2D level. Electrophysiological variability was introduced by using a population of models varying different ionic conductances and extracellular ionic concentrations on a cellular model of chronic AF human atria cardiomyocyte and the anatomical variability was implemented by using models with two different tissue sizes.

This study revealed that simulated dilated atria presented more profiles maintaining reentry, therefore confirming the hypothesis that larger tissues are more prone to fibrillate, an effect that has already been shown at the clinical level (Zou et al., 2005; Qureshi et al., 2014).

The population of models was not only evaluated at basal conditions but also under the effect of three different cardiovascular drugs: flecainide, verapamil, and isoproterenol. Overall pharmacological effects of the implemented drugs matched their clinical characteristics, mainly exhibiting antiarrhythmic effects in the case of flecainide and verapamil (Klein et al., 1979; Echt and Ruskin, 2020). It is interesting to point out that, in the case of dilated tissue experiments, the proportion of non-inducible AF in simulations with flecainide or verapamil was lower than for the smaller tissue size profiles. Furthermore, whereas flecainide and verapamil exhibited mainly an antiarrhythmic effect, isoproterenol exhibited a proarrhythmic effect (Oral et al., 2008). However, some antagonistic effects were observed as exhibited in Figures 3, 4, following the same behavior observed at a clinical level where undesirable effects of antiarrhythmic compounds have been identified in a minority of patients (Nogales Asensio et al., 2007). This confirms our hypothesis that the effect of the drug can be different depending on the specific expression of the ionic currents. Antagonistic effects of antiarrhythmic drugs have previously been studied (Donnelly, 2004) stating the need of identifying the different factors that can lead to this response such as genetics or drug dynamics, which are not considered in this study. However, none of these approaches analyzed the effects on a population of models of the atria.

Furthermore, we observed that the profiles that changed its behavior when adding a drug were similar in the case of the antiarrhythmic drugs simulated in this study. This suggests that specific groups with similar characteristics have analogous responses.

As the number of simulations reached a significant number of samples with an increasing information volume, AI algorithms were applied to reveal patterns in the data. AI application in clinical environments is exponentially increasing and leading toward new diagnostic and treatment techniques (Feeny et al., 2020; Sánchez de la Nava et al., 2021). This trend has also been implemented in the electrophysiology field (Muffoletto et al., 2021; Siontis and Friedman, 2021), in which the use of algorithms has been used for detecting or evaluating proarrhythmicity (Shao et al., 2018; Halfar et al., 2021), classifying different rhythms (Wasserlauf et al., 2019) or automatizing tasks as segmentation (Yang et al., 2017). Moreover, its use in safety pharmacology could be applied to analyze all the data produced by in silico simulations. Particularly, the Random Forest algorithm grouped similar profiles with the same outcome, therefore implementing an AI-driven algorithm able to predict, based on the ionic combinations of each profile, the probability of AF inducibility with excellent predictive values (Sánchez de la Nava et al., 2021).

Data interpretability in the AI field has been demonstrated to be important, especially in the clinical field, where the understanding of the patterns found by algorithms is usually described as a black-box that does not allow to evaluate the biomarker identification and the decision outcome (Nicholson Price, 2018). In this case, the methodology used allowed us to analyze high amounts of data and to understand the clinical implications of the characteristics of each cluster. For example, the first variable that initiated the classification on the Random Forest algorithm was found to be the inward rectifier K⁺ channel that at the clinical level has a crucial role in controlling the frequency and stability of rotors responsible for AF (Pandit et al., 2005). Moreover, up-regulation of I_K1 increased the ability to sustain faster and longer-lasting reentry, predisposing to the development of tachyarrhythmias (Noujaim et al., 2007).

Thus, the identification of this current as the first of the algorithm denotes the importance in the arrhythmia induction and maintenance mechanisms, as overexpression of this repolarization current has been associated with rotor acceleration as a consequence of a reduction in the action potential duration (Atienza et al., 2006; Pandit et al., 2011). The identification of this variable by the algorithm shows that the patterns are explainable and represent clinical scenarios that can be found in a real-life AF population.

This variable was followed by a succession of combinations resulting in several clusters with similar profiles associated with AF inducibility. Interestingly, from the overall implementation of the AI algorithm, two patterns were identified showing an increase in AF inducibility related to changes in extra and intracellular potassium levels and allowed clinical interpretability of the results that facilitate the understanding of the patterns found by the algorithm. A first cluster was found showing increased expression and concentration of potassium ions, which can be directly related to hyperkalemia. Another different pattern observed in a group with different profiles presented low values of both ionic expression and ionic concentration of potassium, which can be directly related to hypokalemic conditions (Guo et al., 2009). At the clinical level, these two conditions have been associated with the induction and triggering of different arrhythmias (Pandit et al., 2011; Skogestad and Aronsen, 2018; Rivera-Juárez et al., 2019; Rakisheva et al., 2020). More in detail, studies have shown that hypokalemia is an independent predictor of developing AF (Krijthe et al., 2013) and that specific cases of AF patients have been identified in whom hyperkalemia could induce malignant arrhythmias (Yan et al., 2018).

From the two main identified pathways that led to sustained AF reentry, the first cluster was associated with increased I_K1, decreased diffusion, decreased extracellular potassium, and overexpression of sodium channels (Figure 5). The identified decreased diffusion can be directly related to lower conduction velocity, which has already been reported as a trigger for re-entrant foci and arrhythmia induction (King et al., 2013). Beyond the key role of I_Na determining excitability (Bezzina et al., 2001), the identification of this first cluster confirms the strong interaction between the molecular correlates of I_Na and I_K1 as part of a common macro-molecular complex, where resting membrane potential hyperpolarization indirectly affect rotor frequency by modifying I_Na availability (Milstein et al., 2012; Ponce-Balbuena et al., 2018).

In the second cluster, both the ionic concentration of potassium and decreased potassium conductivity and extracellular concentration were combined with increased calcium dynamics (Current 1 and Current 2). These currents have a major role during the plateau phase of the action potential and its overexpression contributes to faster repolarization rates, shortening the action potential duration and increasing the risk of early afterdepolarizations which can give rise to the initiation of rotors and fibrillation (Cerrone et al., 2007; Nattel and Dobrev, 2012).

Results presented in this manuscript were based on a population of models of 149 subjects and further samples and variation ranges should be included to explore a wider population. Besides, simulations on the planes do not reflect all proarrhythmic areas present in 3D structures such as pulmonary veins or information including fiber orientation. The use of two different plane sizes did show that the arrhythmia inducibility was higher in bigger tissues but the simplicity of the model presenting constrained borders restricts the overall interpretability of the results, as the real atria contain highly complex structures that play important roles in the initiation and maintenance of AF. In addition, arrhythmia maintenance was considered for planes maintaining rotational activity for more than one cycle, but further analysis should be conducted to evaluate if self-termination occurred in a part of the profiles. Moreover, fibrotic tissue is a relevant condition that predisposes to AF and should be included in future studies and, if possible, the variables that conform to the AI algorithm should be transferable to the clinical practice using metrics such as conduction velocity or rotor dynamic biomarkers. In addition, machine learning algorithms will perform better and will show more robustness with a higher number of samples and different concentrations of each of the drugs. Finally, a greater number of drugs could be simulated in order to enlarge the field of application in which the algorithm can be used.

The identification of specific groups or ionic characteristics that present proarrhythmic effects can be critical in the understanding of new tools for future pharmacological development. At the clinical level, this type of analysis can help to personalize the pharmacological treatments for each of the patients, therefore avoiding possible adverse effects. This study presents, as a result, a new trained algorithm that includes both anatomical and electrophysiological data to evaluate arrhythmia inducibility that is presented as a proof of concept for drug effect evaluation in AF, identifying similar results when compared to the clinic where dilated atria and specific cases of adverse drug effects (Bassett et al., 1997). However, this study includes the implementation of the algorithm based on parameters that are currently difficult to obtain, therefore limiting its current application at the clinical level. New biomarkers describing the ionic characteristics at the patient level should be obtained in order to apply this algorithm.