In Silico Assessment of Class I Antiarrhythmic Drugs Effects on Pitx2-induced Atrial Fibrillation: Insights from Populations of Electrophysiological Models of Human Atrial Cells and Tissues

Electrical remodelling as a result of the homeodomain transcription factor 2 (Pitx2)‐ dependent gene regulation was linked to atrial fibrillation (AF) and AF patients with single nucleotide polymorphisms at chromosome 4q25 responded favorably to Class I antiarrhythmic drugs (AADs). The possible reasons behind this remain elusive. The purpose of this study was to assess the efficacy of AADs disopyramide, quinidine, and propafenone on human atrial arrhythmias mediated by Pitx2-induced remodelling, from a single cell to the tissue level, using drug binding models with multi-channel pharmacology. Experimentally calibrated populations of human atrial action potential (AP) models in both sinus rhythm (SR) and Pitx2-induced AF conditions were constructed by using two distinct models to represent morphological subtypes of AP. Multi-channel pharmacological effects of disopyramide, quinidine, and propafenone were considered. Simulated results showed that Pitx2-induced remodelling increased maximum upstroke velocity (dVdtmax) and conduction velocity (CV), and decreased AP duration (APD) and wavelength (WL). At the concentrations tested in this study, these AADs decreased dVdtmax and CV and prolonged APD in the setting of Pitx2-induced AF. Our findings of alterations in WL indicated that quinidine and disopyramide may be more effective against Pitx2-induced AF by prolonging WL.

A population-based study assessed the influence of AF-associated loci on the response to antiarrhythmic drug (AAD) therapies and showed that carriers of the variant allele at rs10033646 on chromosome 4q25 (Pitx2) responded favorably to Class I AADs [34]. Class I AADs used in AF include flecainide, disopyramide, quinidine and propafenone [35]. The pharmacological effects of flecainide in Pitx2-induced AF were investigated by using a multi-scale computational model and simulated results demonstrated that flecainide is effective for the treatment of Pitx2-induced AF patients by preventing spontaneous calcium release [36]. However, the efficacy of other Class I AADs disopyramide, quinidine, and propafenone on human atrial arrhythmias mediated by Pitx2-induced remodelling remains elusive.
Here, following the Quantitative Systems Pharmacology Framework developed by Ni et al. [55], we constructed populations of in silico models calibrated to values of action potential (AP) biomarkers reported in an experimental dataset on AP recordings [59]. Using these experimentally calibrated models, we simulated and assessed actions of Class I AADs on human atrial electrophysiology. AP duration (APD) and maximum upstroke velocity (dVdtmax) were quantified to evaluate anti-AF effects of Class I AADs on AP at the cellular level, while conduction velocity (CV) and wavelength (WL) were quantified to assess the effects of Class I AADs on AP propagation at the tissue level. Sensitivity analyses of AP biomarkers were applied to understand the ionic mechanisms underlying Pitx2-induced AF and the efficacy of these AADs. Finally, we performed populationbased simulations of Pitx2-induced remodelling and predicted a reduction in APD and WL and an increase in CV and dVdtmax. Further simulations of actions of Class I AADs on Pitx2-induced AF exhibited APD prolongation and a reduction in CV and dVdtmax. Our results showed that quinidine and disopyramide led to WL prolongation compared to the drug-free AF conditions, but propafenone caused WL shortening. These findings suggest that quinidine and disopyramide may be effective against Pitx2-induced AF. Figure 1. Flow chart illustrating the process for in silico assessment of Class I antiarrhythmic drugs effects on Pitx2-induced atrial fibrillation (AF). Note: There are two distinct morphological subtypes of human atrial action potential (AP) [59] and thereby the Bai et al. model with a notch-and-dome AP morphology [27] and the Grandi et al. model with a triangular AP shape [60] were used in the present study. Class I antiarrhythmic drugs assessed include quinidine, disopyramide and propafenone. Abbreviations: Pitx2-The homeodomain transcription factor 2, SR-Sinus rhythm; dVdtmax-Maximum upstroke velocity, RMP-Resting membrane potential, APD50 and APD90-AP duration at 50% and 90%, respectively, CV-Conduction velocity and WL-Wavelength. Table 1. Parameters associated with ionic properties for constructing populations of atrial models.
A flow chart illustrating the process for in silico assessment of Class I AADs effects on Pitx2induced AF is presented in Figure 1. First, the initial population of sinus rhythm (SR) models was created by randomly perturbing parameters associated with ionic properties ( Table 1) of the baseline human atrial cell models [41]. The Bai et al. model displaying a type-1 AP with notch-and-dome morphology [27] and the Grandi et al. model displaying a type-3 AP with typical triangular shape [60] were chosen as the baseline human atrial cell models to represent distinct morphological subtypes of human atrial AP [62]. Parameters associated with ionic properties were allowed to vary independently according to a log-normal distribution and sigma was set to be 0.2 to cover a range of variability similar to that seen in experiments based on previous studies [41,55,63,64]. Second, we used the initial models to simulate human atrial AP by considering stimulation frequency (1Hz) under the experimental conditions [59] and calculate AP biomarkers of each initial model. The initial models generated in the previous step were selected to constitute the SR population whose simulated electrophysiological properties are in range with the same properties in experimental data on AP biomarkers (including dVdtmax, RMP, APD50 and APD90) in Table 2. This step yields the experimentally-calibrated population of SR models [54]. Third, electrical remodelling due to impair Pitx2 ( Table 3) was introduced into SR model variants to generate the initial population of AF models. Fourth, we used experimental AP biomarkers to calibrate the AF population by excluding the model variants in which values of simulated biomarkers were outside the experimentally observed range ( Table 3) reported by Ravens et al [59]. Fifth, the blocking effects of AADs on ion channels ( Table 4) were incorporated into experimentally-calibrated Pitx2-induced AF models to evaluate their effects on the virtual atrial myocytes. Finally, sensitivity analyses [41,55] of dVdtmax and APD90 at the cellular level and CV and WL at the tissue level were applied to understand modulations by ionic parameters of Pitx2-induced remdoelling and ion channels affected by AADs. On the basis of model responses, AADs that eliminates the arrhythmogenic propensity of the atrial substrate arising from Pitx2induced remodelling were selected. Note: Antiarrhythmic drugs investigated include disopyramide, quinidine and propafenone. Ion currents associated with these drugs included INa, ICaL, Ito, IKs, IKr, IKur and IK1. Blocking effects of drugs on ion currents were modelled with the half-maximal inhibitory concentration (IC50) and Hillcoefficient (nH) value. IC50 and nH were extracted from the literature. Note: Values are presented in mean ± standard error. AP biomarkers include APD50, APD90, RMP and dVdtmax. Abbreviations: dVdtmax-Maximum upstroke velocity, RMP-Resting membrane potential, and APD50 and APD90-AP duration at 50% and 90%, respectively.

The experimentally-calibrated populations of human atrial myocytes under SR and AF conditions
Representing a type-1 AP with notch-and-dome morphology, the Bai et al. model [27] was used to generate the initial SR population of human atrial cell model variants. According to the physiological range of AP biomarkers (dVdtmax, RMP, APD50 and APD90 values at 1Hz pacing) measured experimentally ( Table 2), populations of 610 SR models out of the initial pool of 1,200 models were selected to generate the experimentally-calibrated population of SR models. Next, we incorporated Pitx2-induced electrical remodelling ( Table 3) into the SR population of 610 variants to generate the initial AF population, and then calibrated this AF population to the experimentally measured AP biomarkers ( Table 2). Representative AP traces for the experimentally-calibrated SR population of 610 models and AF population of 214 models, respectively, are shown in Figure 2a and Figure 2b. Distributions of AP biomarkers of the experimentally-calibrated SR models are compared to that of AF models (Figure 2c-d). In detail (  Representing a type-3 AP with a typical triangular shape, the Grandi et al. model [60] was also used to generate the experimentally-calibrated populations according to the physiological range of AP biomarkers [59] measured experimentally ( Table 2). 153 SR models and 79 AF models were selected to constitute the SR population and the AF population, respectively. Figure 3 shows representative AP traces (Figure 3a

Sensitivity analysis revealed alterations in depolarization and repolarization of AP
Experimental studies have demonstrated that sensitivity analyses of dVdtmax and APD90 can provide insights into alterations in the depolarization rate [78][79][80][81] and the repolarization time [82] of the cardiac tissue. We applied partial correlation analysis to predict alterations in dVdtmax and APD90 when electrical remodelling (listed in Table 3) due to impaired Pitx2 was introduced or ion channels (listed in Table 4) were blocked by Class I AADs. Parameters associated with Pitx2-induced electrical remodelling include GNa, GCaL, GKs, GK1, Grel and Gup, while parameters affected by actions of AADs (including disopyramide, quinidine and propafenone) include GNa, GCaL, Gto, GKs, GKr, GKur and GK1. We used partial correlation analysis to calculate partial correlation coefficients (PCCs) to quantify correlations between the parameter values and dVdtmax and APD90, respectively. In  Table 1) are shown in Figure  1S and Figure

Antiarrhythmic effects of Class I drugs on dVdtmax and APD90 at the cellular level
According to experimental data on actions of Class I AADs on ion currents, their effects were incorporated into experimentally-calibrated models in the AF population created with the Bai et al. model. Parameters associated with the effects of AADs were GNa, GCaL, Gto, GKs, GKr, GKur and GK1 ( Table  4). Class I AADs investigated here included disopyramide, quinidine and propafenone. Taking into account plasma protein binding, estimates of the most likely unbound concentrations of propafenone, disopyramide and quinidine have been given as ~0. 15 [37,87]. Figure 5 shows the actions of disopyramide, quinidine and propafenone on dVdtmax and APD90 in the AF condition. It can be seen in Figure 5a that disopyramide, quinidine and propafenone reduced dVdtmax in a dose-dependent manner, with propafenone decreasing dVdtmax to a greater extent than disopyramide and quinidine. Values of dVdtmax upon application of disopyramide and quinidine were smaller than those in the drug-free AF condition, but were larger than that in the drug-free SR condition. However, values of dVdtmax upon application of propafenone were smaller than those under drug-free AF and SR conditions. All drugs also prolonged APD90 in a dose-dependent manner (Figure 5b), with quinidine producing a slightly larger increase in APD90 across all concentrations investigated. In detail (  propafenone decreased dVdtmax to a greater extent than disopyramide and quinidine. Values of dVdtmax upon application of propafenone and quinidine in the AF condition were smaller than that in the drug-free SR condition, but dVdtmax upon application of disopyramide in the AF condition were larger than that in the drug-free SR condition (Figure 6a). It can also be seen in Figure 6b that APD90 upon application of disopyramide in the AF condition were larger than that in the drug-free SR condition, whereas values of APD90 upon application of propafenone and quinidine in the AF condition were smaller than that in the drug-free SR condition.

Antiarrhythmic effects of Class I drugs on CV and WL at the tissue level
To evaluate the effects of Class I drugs on CV and WL, we constructed one-dimensional (1D) models of human atrial strands with the Bai et al. model to investigate the responses of electrical waves to AADs in tissue. Figure 7 shows CV and WL upon application of disopyramide, quinidine and propafenone in the AF condition, compared with those in the drug-free SR and AF conditions. CV was increased from 0.56 ± 0.09 mm/ms for the SR condition to 0.58 ± 0.09 mm/ms for the AF condition. It can be seen in Figure 7a that disopyramide, quinidine and propafenone reduced CV in a dose-dependent manner, with quinidine and propafenone decreasing CV to a greater extent than disopyramide. Values of CV upon application of disopyramide and quinidine were smaller than those in the drug-free AF condition, but were larger than those in the drug-free SR condition. However, values of CV upon application of propafenone were smaller than those under drug-free AF and SR conditions ( Table 6). Figure 7b shows values of WL under drug-free SR and AF conditions, and upon application of various concentrations of disopyramide, quinidine, and propafenone. WL was decreased from 141.05 ± 27.93 mm for the drug-free SR condition to 108.15 ± 20.19 mm for the drug-free AF condition. Compared with the drug-free AF condition, propafenone shortened WL in a dose-dependent manner, whereas disopyramide and quinidine prolonged WL in a dose-dependent manner. Quinidine prolonged WL to a greater extent than disopyramide. And upon application of 5μM quinidine (Quin_H) in the AF condition is close to that in the drug-free SR condition (132.22 ± 24.13 mm vs. 141.05 ±27.93 mm). In order of effects of AADs on WL, these drugs are quinidine, disopyramide and propafenone.  (Figure 8a). In addition, values of CV upon application of various concentrations of propafenone were smaller than those under drug-free SR conditions. Figure 8b shows that disopyramide prolonged WL in a dosedependent manner, whereas quinidine and propafenone shortened WL in a dose-dependent manner. However, values of WL upon application of quinidine (Quin_L and Quin_M) and propafenone (Prop_L) were larger than those under the drug-free AF condition, indicating that quinidine is more effective that propafenone. In order of effects of AADs on WL, these drugs are disopyramide, quinidine and propafenone.
Quantitative summaries of the effects of Class AADs on human atrial electrical activity in the Pitx2-induced AF condition are listed in Table 6 for the Bai et al. model and Table 7 for the Grandi et al. model.

Discussion
Population-based studies have assessed the influence of common single-nucleotide polymorphisms related to AF on the response to AAD therapies and showed that carriers of the variant allele at rs10033646 on chromosome 4q25 responded favorably to the Class I AADs [34]. The actions of Class I AADs disopyramide, quinidine, and propafenone were assessed in the context of Pitx2-induced AF using a population-based Quantitative Systems Pharmacology Framework [55]. Through sensitivity and statistical analyses of our atrial cell and tissue populations, we found that WL prolongation could be achieved by quinidine and disopyramide. This study provides clinicallyrelevant insights into the pharmacology of WL prolongation by evaluating and comparing the actions of all three drugs in the context of Pitx2-induced AF, offering an important step toward in silico optimization of pharmacological therapy in this context.

Main findings
The major findings presented in this study are as follows. (1) Populations of models based on two human atrial AP models are able to mimic a wide range of inter-subject variability in human atrial AP properties as exhibited in a set of AP measurements from over 379 SR and AF patients. (2) Pitx2-induced remodelling (as reported in [36]) predicts abbreviated APD90 and increased dVdtmax at the cellular level, and increased CV and shortened WL at the tissue level in SR versus AF conditions, as reported in experimental studies [17,59] (Table 2) AADs disopyramide, quinidine, and propafenone prolonged APD90 and decreased dVdtmax and CV in a dose-dependent manner in drug-free versus drug-bound AF conditions. (5) Disopyramide and quinidine prolonged WL, whereas propafenone shortened WL. Disopyramide and quinidine are more effective against Pitx2-induced AF than propafenone. A summary of findings regarding proarrhythmic mechanisms of Pitx2-induced AF and anti-arrhythmic actions of selected Class I drugs is given in Figure 9.

Changes in AP are linked to Pitx2-induced remodelling
Research studies on Pitx2-induced AF have demonstrated that Pitx2-induced remodelling contributes to atrial cellular electrophysiological changes in AF patients [11,17,26]. Experimental data showed that alterations in human atrial AP included AP shortening and increased upstroke velocity [59] and our simulated effects of Pitx2-induced remodelling on AP ( Table 2 and Table 5) are concordant with these experimental findings.
Previous experiments on isolated human atrial myocytes have demonstrated that Pitx2-induced remodelling of ion channels, particularly for IKs and ICaL, may contribute to the clinically significant association between impaired Pitx2 and AF [26]. Further simulations in our previous study indicated the Pitx2-induced changes in IKs and ICaL led to APD shortening, facilitating sustained re-entry in 3D anatomical atrial geometry [30]. Our simulations in the Pitx2-induced AF population of models show that AP shortening is determined by ICaL and IKs, and APD is negatively correlated with IKs and positively correlated with ICaL (Figure 4). This is consistent with findings from previous studies that upregulation of IKs and downregulation of ICaL due to Pitx2-induced remodelling critically contribute to the abbreviation of APD [30,33].
Furthermore, increased upstroke velocity in Pitx2-induced AF may be resulted from remodelled INa arising from impaired Pitx2. Previous experiments in atrial cardiomyocytes of Pitx2-deleted mice showed increased action potential amplitude and shortened APD [17]. Further gene analysis showed that reduced Pitx2 dose caused an increase in the expression of SCN5A [17] encoding the alpha subunit of the human cardiac voltage-gated sodium channel (INa). In addition, an increase in the mRNA level of SCN1B encoding a beta-1 subunit of voltage-gated sodium channels (INa) was observed in patients with familial AF arising from the Pitx2c mutation p.Met207Val [61]. Decreased mRNA levels of SCN5A and SCN1B were also observed in atrial-specific Pitx2 mutant mice [22,25], but changes in T-box transcription factor Tbx5 was not investigated and downregulation of sodium channel genes may be resulted from impaired Tbx5. Previous experiments in atrial cardiomyocytes of Tbx5-deleted mice showed reduced Pitx2 and sodium channel genes although reduced Pitx2 without alterations in Tbx5 caused an increase in sodium channel genes [17]. Based on these data on remodelled INa due to impaired Pitx2, Pitx2-induced AF model was developed and predicted an increase in action potential amplitude and maximum upstroke velocity [36]. Consistent with these findings, our simulations in the Pitx2-induced AF population of models show that increased upstroke velocity is determined by INa and is positively correlated with INa (Figure 4). Further simulations at the tissue level predicted that increased INa led to fast conduction of AP propagation.

Anti-arrhythmic effects of Class I drugs in Pitx2-induced AF
Despite the prevalence of AF and decades of research, antiarrhythmic therapies for AF continue to have limited efficacy and safety [35]. A population-based study assessed the influence of AFassociated loci on the response to antiarrhythmic drug therapies and showed that carriers of the variant allele at rs10033646 on chromosome 4q25 (Pitx2) responded favorably to Class I AADs [34]. And Class I AADs frequently used clinically include flecainide, disopyramide, quinidine and propafenone [35]. Anti-arrhythmic effects of flecainide on Pitx2-induced AF were investigated in our previous study and simulation results indicated that flecainide has antiarrhythmic effects on AF due to impaired Pitx2 by preventing spontaneous calcium release and increasing wavelength [36]. In the present study, the effectiveness of disopyramide, quinidine and propafenone was assessed with WL. Simulated results showed that disopyramide and quinidine were moderately effective by prolonging WL, whereas propafenone was shown to be ineffective by abbreviating WL in Pitx2-induced AF ( Table 6). WL is the product of APD90 and CV. Increased APD90 ( Figure 5) and decreased CV (Figure  7) occurred upon the application of disopyramide, quinidine and propafenone. Under the application of disopyramide or quinidine, the extent of AP prolongation is larger than the degree of CV reduction, resulting in an increase in WL. In contrast, the extent of CV reduction is much larger than the degree of AP prolongation, leading to a reduction in WL. Alterations in WL can be linked to block effects of Class I drugs on ion channels ( Table 4) which included GNa, GCaL, Gto, GKs, GKr, GKur and GK1. Under the application of Class I drugs in Pitx2-induced AF, inhibition of sodium currents would decrease dVdtmax and thereby reduce CV, while block of potassium currents would prolong APD90. Therefore, the antiarrhythmic effects of Class I drugs in Pitx2-induced AF can be attributed to their potent block of sodium and potassium channels. And the effectiveness of drugs is determined by the extent of their block effects on sodium and potassium channels. Therefore, disopyramide and quinidine are more effective than propafenone in Pitx2-induced AF and this is because the extent of their block effects on potassium channels is larger than that on sodium channels. Collectively, our results demonstrate that a combined block of sodium and potassium currents can exert synergistic antiarrhythmic effects and therefore is a valuable therapeutic for Pitx2-induced AF.

Model-and AP-shape independence of prediction of the antiarrhythmic effects of Class I drugs
It is well known that human atrial myocytes display two distinct AP morphologies: the type-1 AP shows a notch-and-dome morphology and the type-3 AP shows a typical triangular morphology [59,62]. Therefore, Bai et al. model displaying the type-1 AP [27] and Grandi et al. model displaying the type-3 AP [60] were chosen to generate populations of human atrial AP. According to the stimulation protocol (1Hz) used in experiments [59], we ran simulations and reproduced Pitx2induced AP morphology, with both dVdtmax increase and AP shortening being reported experimentally [59]. Further simulations predicted the antiarrhythmic effects of Class I drugs on APD90 and dVdtmax, but their effects are dependent on baseline AP morphology [88,89] seen in previous modelling [55]. While some of the differences in model responses might be related to distinct AP morphologies, we cannot exclude model-dependencies due to distinct cellular model structure and models of ionic and calcium handling processes. Interestingly, analyses using both models demonstrated that these Class I drugs disopyramide, quinidine and propafenone prolonged APD90 and reduced dVdtmax in a dose-dependent manner. Therefore, these results demonstrated that these Class I drugs consistently produced anti-AF effects independent of the baseline electrophysiological characteristics.
Although Class I drugs disopyramide, quinidine and propafenone produced similar anti-AF effects (AP prolongation and dVdtmax and CV reduction), the degree of changes in APD90 and dVdtmax is different between Bai  Different alterations in CV and APD90 leading to different changes in WL, which can be used as one of the indexes for evaluating the antiarrhythmic effects of drugs at the tissue level [55]. If the extent of AP prolongation is larger than that of CV reduction, WL is prolonged. Conversely, WL is shortened.

Limitations and Future Work
Several limitations specific to this study are addressed here. Firstly, the electrophysiological representation of AF-induced remodelling in the human atrial AP model is based on data from previous mice models of Pitx2-induced AF, however, because of the lack of experimental data on humans. Special attention must be paid to the differences between mice and human atrial cells [90]. Further, we simulated AF model populations with all identified targets associated with impaired Pitx2 [11,17,18,22,23,26]. However, Pitx2-induced remodelling is different in published experimental studies [11,17,18,22,23,26], atrial cell types [11], AF stages [91] and AF patients with different age [92]. The complexity of Pitx2-induced remodelling may present complex responses to these Class I drugs and require more realistic heterogenous descriptions in cellular and tissue simulations. Secondly, effects of these Class I AADs on ion currents were modelled by changing maximum current conductances with a simple pore block scheme based on IC50 and nH values, but it might be important to incorporate state-dependence and use-dependence of antiarrhythmic drugs in evaluating realistic compounds [37,63,87]. Model variations for these characteristics (including the affinity of the drug compound to various gating modalities, binding kinetics, drug polarity, charge for the drug-channel interaction and so on) are far beyond the scope of the current study and should be explored in future studies. In addition, IC50 values are chosen based on experimental data from atrial cells (where data are available) and large ventricular cells (where atrial data are not available), but the large variability in IC50 was influenced by various experimental conditions, including different species, different cell types [93], temperature [94], ionic concentrations and voltage-clamp protocols. Therefore, special attention should be paid to explain our simulated results. Thirdly, considering the stimulation protocol used under the experimental conditions [59], we performed all simulations only at 1Hz and rate-dependent APD modulation will be further investigated when experimental data at different frequencies are available. Here, we only used APD50, APD90, RMP and dVdtmax as biomarkers to calibrate populations of human atrial AP models. And parameter unidentifiability, which can be attenuated if broader experimental datasets are available, is a potential limitation. Finally, we analyzed the effects of Class I drugs by quantifying changes in APD90 and dVdtmax at the cellular level and CV and WL in tissue. This approach is consistent with the accepted mechanisms of action of Class I drugs and is limited in that it didn't account for composite metrics, including calcium-associated biomarkers, afterdepolarizations, alternans and AP propagation dynamics. Nevertheless, mechanisms underlying actions of Class I drugs are highly complex and future investigations should be carried out.

Experimantal dataset
The Experimental dataset on AP recordings was used in the present study for calibration of the human atrial electrophysiology. The data set published in the previous study [59] comprised of 480 instances from 379 patients: AP recordings from n = 256 right atrial appendages of N = 221 SR patients and from n = 224 right atrial appendages of N = 158 AF patients. Human myocytes were isolated enzymatically from atrial appendages and APs were recorded with standard intracellular microelectrodes in atrial trabeculae. Preparations were electrically stimulated at a single constant rate of 1 Hz for 60 min with isolated square-wave stimuli of 1 ms duration, two times threshold intensity. Obtained human atrial APs displayed spike-and-dome and a more triangular conformation. The following parameters were quantified to characterize intersubject variability in human atrial AP: dVdtmax, RMP, and APD50 and APD90. Compared with APs of SR myocytes, APD50 and APD90 were reduced and dVdtmax was increased for APs of AF atrial cells ( Table 2). More information (including ethics approval, informed consent and basic information of participants) regarding the experimental conditions under which the data were collected is available in the study of Ravens et al [59].

Mathematical models representing distinct morphological subtypes of human atrial AP
To represent spike-and-dome and a triangular conformation of human atrial APs observed in experiments [59,62], the Bai et al. model displaying a type-1 AP with notch-and-dome morphology and the Grandi et al. model displaying a type-3 AP with typical triangular shape were chosen as a base to construct the computational AP model populations.
The Bai et al. model developed by our group is able to reproduce human AP morphology, APD rate dependence and triggered activity, i.e., early afterdepolarizations (EADs), delayed afterdepolarizations (DADs) and spontaneous depolarizations [27]. This biophysically-detailed model of human atrial cellular electrophysiology also has been used to investigate mechanisms underlying Pitx2-induced AF and here we provide a brief description. It includes representation of the 13 transmembrane ionic currents and 2 main intracellular calcium flows, including fast sodium current (INa), transient outward potassium current (Ito), rapid delayed rectifier potassium current (IKr), slow delayed rectifier potassium current (IKs), ultrarapid delayed rectifier potassium current (IKur), inward rectifier potassium current (IK1), L-type calcium current (ICaL), background sodium current (IBNa), background calcium current (IBCa), plateau potassium current (IPK), plateau calcium current (IPCa), sodium-potassium pump current (INaK), sodium-calcium exchange current (INcx), the calcium flow (Iup) through the sarcoplasmic reticulum calcium ATPase (SERCA) and calcium release flow (Irel).
Comparative simulations were carried out using the human atrial cell model of Grandi et al. [60,95] where (14.838 nS/pF) is the maximal conductance, , ℎ and are three gate variables for INa. is the membrane potential, and is the sodium equilibrium potential. ∞ , ℎ ∞ and ∞ denote steady-state activation, steady-state inactivation and steady-state inactivation, respectively. , ℎ and are the time constants for ∞ , ℎ ∞ and ∞ , respectively. Based on these models of ion channels, the electrophysiological behavior of human atrial cells can be modelled with the following differential equation: where is the membrane potential, is time, is the total sum of transmembrane ionic currents (listed in Table 1), is the externally applied stimulus current with amplitude of -80 / and a duration of 0.5 ms, and is the membrane capacitance (1 ). In order to investigate the electrophysiological behavior of human atrial tissues, a multicellular model of homogeneous atrial tissue with 200 nodes spaced evenly by 0.15 mm was constructed to quantify APD90 and CV for calculating WL. The propagation of APs was governed by the partial differential equation: where is a scalar coefficient describing the intercellular electrical coupling via gap junctions. In simulations, was set to be a constant value of 0.154

Modelling Pitx2-induced AF
To obtain human AF myocytes models that reproduced the experimentally observed changes in the mRNA levels corresponding to key proteins under Pitx2-induced electrical remodelling conditions, we assumed that these changes in mRNA expression are quantitatively reflected at the final functional level of ion channels [29] and incorporated alterations in the maximal conductances of ionic currents due to Pitx2-induced electrical remodelling into the Bai et al. (or the Grandi et al.) model. These Pitx2-induced changes in the ionic channel properties have been well characterized in many experimental studies [17,18,[23][24][25][26]61]. However, the identified targets were different among these studies: while several studies showed no Pitx2-associated changes [18] and reductions [22,25] in sodium channel gene expression, significant increases [17,61] in this channel gene expression have also been well-documented. Similarly, calcium channel gene expression was found decreased in this study [12], but increased in others [22,23,25,26]. Further, whereas overexpression of potassium gene encoding IK1 was identified in the study of Kirchhof et al. [11], underexpression of this gene was observed in some studies [18,22]. In addition, all these studies have shown upregulation of genes encoding IKs [11,12,15,26,61], Irel [12,17,23,25]and Iup [22,24,67,70,71]. Modifications made for each model are summarized in Table 3.

Simulations of actions of Calss I AADs in Pitx2-induced AF
In our previous studies [37,87], the actions of Class I AADs disopyramide, quinidine and propafenone on human atrial electrophysiology were simulated in the setting of hERG-linked short QT syndrome [37]. The effects of these Class I AADs on ion currents were modelled by changing maximum current conductances with a simple pore block scheme based on IC50 and nH [31]. Maximum current conductances associated with these drugs included GNa, GCaL [83][84][85].

Generation and calibration of populations of models
To capture intersubject variability, populations of sampled models of human atrial electrophysiology for SR were generated based on the Bai et al. and the Grandi et al. models. All models in each population shared the same equations but parameters of ionic current conductances in determining the human atrial AP were varied with respect to their original values. These parameters were independently varied following a log-normal distribution and sigma was set to be 0.2 to cover a range of variability similar to that seen in experiments based on previous studies [41,55,63,64]. The size of the SR population was set to be 1200 for convergence of the sensitivity coefficients [101]. Following the ASME V&V40 Standard proposed by the Subcommittee of the American Society of Mechanical Engineers (ASME) on Verification and Validation (V&V) in Computational Modeling of Medical Devices [102] for developing a structured approach for establishing the credibility of computational models for a specific use, we used the candidate models to simulate human atrial AP by considering stimulation frequency (1Hz) under the experimental conditions [59] and calculate AP biomarkers of each candidate model. The candidate models generated in the previous step were selected to constitute the SR population whose simulated electrophysiological properties are within 3 standard error from the mean experimental values in SR patient data on AP biomarkers (including dVdtmax, RMP, APD50 and APD90) in Table 2. This step yields the experimentally-calibrated population of models. Modifications to ionic currents due to Pitx2induced remodelling ( Table 3) were introduced into SR model variants to generate the candidate AF models. These candidate AF models were calibrated to 3 standard error from the mean experimental values in AF patient data ( Table 2). Then, actions of Class I AADs on ion channels were incorporated into the experimentally-calibrated AF models to investigate their effects on dVdtmax, APD90, CV and WL.

Simulation protocols
The stimulation protocol mimics the one used by Ravens et al. [59] to obtain AP measurements in SR and AF cardiomyocytes. Cell mathematical models were initially preconditioned by pacing at a basic cycle length of 1000 ms until the steady-state was reached. A stimulus with an amplitude of −80 / and a duration of 0.5 ms is applied at each pace. Considering computational cost for reaching steady-state, a series of 10 conditioning waves was initiated by supra-threshold stimuli with an amplitude of −80 / and a duration of 1.0 ms to the 3 nodes at the strand end. The last beat was recorded for APD90 and CV analysis.

Sensitivity analysis
To quantify the relative importance of ionic conductances in determining changes in AP biomarkers, PCCs were used on the SR and AF populations to evaluate the role of each ionic current [101]. Partial correlation is a method to find correlations between two variables, after accounting for the linear effects of one or more additional variables [103]. PCC between and , given the set of additional variables , is then defined as the correlation coefficient between the residuals = −̂ and = −̂ [101]. ̂ and ̂ are the respective sample means or the following linear regression models: where ( , ) represents the covariance between and , while ( ) and ( ) are respectively the variance of and variance of .

Software, numerical methods and statistical analysis
The Bai et al. (available from the repository CellML http://models.cellml.org/workspace/520)and the Grandi et al. models (freely available at https://github.com/drgrandilab) was implemented in MATLAB 2018a (The MathWorks, Natick, MA, USA) using the stiff ordinary differential equation solver ode15s and analysis of biomarkers was also performed using MATLAB. In addition, the 1D models were implemented in C++. The ODEs were solved using the Forward Euler method. Our user project containing newly created datasets and the simulation codes used in this study is available to download from the GitHub website (https://github.com/aspirerabbit). All simulations and data analyses were performed on a computing cluster with Intel(R) Xeon(R) CPU E5-2690 v4 @ 2.60GHz 32 CPUs 28 CPUs (56 threads) + 128GB. Statistical significance in differences on ionic conductance distributions between populations was evaluated by using the Mann-Whitney U test. A probability <0.05 was considered statistically significant.

Conclusions
In conclusion, populations of models reproduce the variability in human atrial AP properties measured in samples obtained from patients and AF models predict AP shortening and fast conduction in Pitx2-induced remodelling conditions observed in experiments. State-of-the-art Quantitative Systems Pharmacological simulations demonstrated that disopyramide, quinidine and propafenone produce AP prolongation and slow conduction in the setting of Pitx2-induced AF. However, disopyramide and quinidine were more effective in prolonging WL than propafenone.

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