RESEARCH

,


INTRODUCTION
Investigations of cardiomyocyte electrical properties in animal models can be controlled to limit variability; in contrast, heterogeneity in the human population and representativeness of the myocytes obtained from cardiac biopsies are key challenges in our understanding of human cardiac physiology. Ion current density can also vary in response to intracellular and external stimuli. In particular, ionic current properties are modulated by signaling molecules such as nitric oxide (13,14,36), hormones (1,2,47), nutrients (57), the circadian rhythm (19), and temperature (10,11). Therefore, when investigating cardiac physiology, we are facing a moving target that is modulated by a host of internal and external factors.
A further constraint in the characterization of human physiology is limited access to human tissue. For this reason, data sets in cardiac electrophysiology often consist of standard action potential (AP) measurements at a single frequency; sometimes only biomarker values, such as the AP duration (APD), are available (18,20,35,39). Similarly, ionic currents are frequently characterized using voltage-clamp recordings in cells that are different from those used for AP measurements, reflecting poor survival of cells subjected to different experimental configurations, solutions, and pharmacological action. Experimental electrophysiological recordings therefore offer a valuable but restricted snapshot of specific properties in isolated tissue/cells. Importantly, integration of these resources remains an open challenge (38).
Construction of biophysically detailed models and simulations can help integrating experimental data and improving knowledge of human electrophysiology (5,27). Recently, consideration of biological variability has triggered further investigations into its underlying ionic basis and implications for disease and treatment (25,29,40,54). However, previous studies of variability in human atrial electrophysiology only had access to AP biomarkers, measured in heterogeneous patient cohorts (20,39).
Here, we investigated how variability in specific ionic currents impacts the phenotypic variability observed in human atrial APs and Ca 2ϩ transients using three different and wellestablished models of human atrial electrophysiology. To do so, we used a rich ex vivo data set consisting of APs recorded at five pacing frequencies in human right atrial myocytes in addition to measurements of four key currents underlying atrial plateau and repolarization phases [namely, the transient outward K ϩ current (I to ), inward rectifier K ϩ current (I K1 ), L-type Ca 2ϩ current (I CaL ), and atrial-specific ultrarapid K ϩ current (I Kur )] in a controlled cohort of patients under sinus rhythm.
Our results support a wide variability in ionic current densities for both K ϩ and Ca 2ϩ currents in human atrial cardiomyocytes. In addition to the role of I CaL in modulating total intracellular Ca 2ϩ , our findings suggest that differences in human atrial Ca 2ϩ transients are primarily driven by intracellular Ca 2ϩ -handling and Ca 2ϩ membrane extrusion. By exploiting the constructed populations of human atrial models, we demonstrate that diastolic Ca 2ϩ fluctuations are determined not only by ryanodine receptor (RyR) or sarco(endo)plasmic reticulum Ca 2ϩ -ATPase (SERCA), as previously reported (33), but importantly also by the Na ϩ /Ca 2ϩ exchanger (NCX). By integrating experiments and simulations, we extend the existing knowledge on ionic and subcellular mechanisms underpinning variability in human atrial electrophysiology, with important implications for further understanding the mechanisms of intracellular atrial Ca 2ϩ dysregulation.

METHODS
Experimental data set. Whole cell current-and voltage-clamp experiments were used to measure APs and I to, IKur, ICaL, and IK1 densities in isolated human atrial myocytes. Myocytes were obtained from the right atrial appendage of patients in sinus rhythm, aged 69 Ϯ 10 yr and undergoing coronary artery bypass graft or aortic valve replacement (Table 1). APs were measured in n ϭ 29 cells at five pacing frequencies (0.25, 0.5, 1, 2, and 3 Hz). Figure 1A shows the AP biomarkers computed in this study, that is, APDs at 20%, 50%, and 90% repolarization (APD 20, APD50, APD90, respectively), resting membrane potential (RMP), and AP amplitude (APA). In voltage-clamp experiments, pharmacological inhibition and/or voltage steps were used to isolate the desired components of the total transmembrane current. IKur, Ito, ICaL, and IK1 densities were measured in n ϭ 30, 30, 18, and 35 cells, respectively. Figure 1, B-E, shows our ex vivo data set at the level of the AP and current densities. A detailed description of experimental protocols is provided in Experimental Data Set in the APPENDIX. Investigations were approved by the Research Ethics Committee (reference number: 07/Q1607/38), and investigations were conducted in accordance with the principles of the Declaration of Helsinki. All patients gave informed written consent.
Constructing in silico populations of human atrial electrophysiological models. Ex vivo variability in AP and ionic current data cannot be captured with a single in silico model of human atrial electrophysiology ( Fig. 1, B-E), motivating the use of a population-of-models approach (4). For this purpose, we select the following three generic in silico models from the literature based on the degree of similarity between model outputs and our ex vivo data (Fig. 1, B-C): 1) the Maleckar et al. model (23) yields a triangular AP with APD closest to our experimental data; 2) the Grandi et al. model (15) produces a triangular AP; however, its APD is longer than in the Maleckar et al. model and lies further away from experiment; and 3) the AP of the Courtemanche et al. model (7) is characterized by a spike-and-dome morphology with the longest APD out of the three models and is the furthest in AP shape and biomarkers from our ex vivo data.
For simplicity, the models are henceforth referred to by their first authors' names. Interestingly, despite substantial differences in the three models' AP properties, the ionic currents thought to be the major determinants of the human atrial AP, i.e., Ito, IK1, ICaL, and atrialspecific I Kur, are not only similar between the models but are also within/close to the experimentally observed ranges of variability in our ex vivo data set (Fig. 1, D and E).
In our experimental data set, peak ICaL density occurred at 10 or 20 mV in n ϭ 10 and 8 cells, respectively (Fig. 1E). Although positively shifted, this is in agreement with reported variability in voltages of peak ICaL in human atrial myocytes (28,30,37). To capture this in silico, modifications in I CaL gating were introduced in the three generic models. Details are provided in Generic Models of Human Atrial Electrophysiology in the APPENDIX.
To capture cell-to-cell variability as seen in our experimental data set, we constructed populations of in silico models sharing the same equations but with different parameter sets and in range with experimentally observed variability (4). For every generic model, a corresponding population of models was constructed. As previously discussed (for a review, see Ref. 25), our underlying hypothesis was that variability in ion channel densities is the main determinant of AP variability in cellular electrophysiology. In all three generic in silico models, maximal conductances/permeabilities of Ito, IKur, ICaL, rapid (I Kr) and slow (IKs) components of the delayed rectifier K ϩ current, I K1, fast Na ϩ current (INa), Na ϩ -K ϩ pump current (INaK), NCX current (I NCX), uptake current (Jup), and release current (Jrel) were varied simultaneously over Ϯ100% of their baseline model values. This choice of currents was based on their influence on human atrial AP morphology and Ca 2ϩ handling following previous investigations and sensitivity analyses (20,39). A wide range of parameter values was used to explore a variety of possible ionic scenarios reflecting the heterogeneity in the human population and our ex vivo data (Fig. 1, D  and E). Such values are supported by physiological evidence in both inward and outward atrial currents in human patients in sinus rhythm (8,48,49). Because there is no experimental evidence on covariation between maximal conductances of two or more ionic currents (12), model parameters were varied independently; should experimental evidence on parameter covariation become available, it can be incorporated into the methodology. To account for heterogeneity in ex vivo peak ICaL density (Fig. 1E), two versions of every generic model, corresponding to the distinct ICaL peaks, were used to build the corresponding in silico population. To verify that the resultant in silico populations were independent from sampling methodology, the following two methods for generating parameter sets from a high-dimensional parameter space were used: Latin hypercube sampling (LHS; see Ref. 24) and sequential Monte Carlo (SMC; see Ref. 9). Under LHS, the range of values that every varied parameter can take is subdivided into N equally spaced intervals, and N parameter sets are generated in such a way that every varied parameter takes a value within the specific interval only one time. We use LHS to create a candidate population of N ϭ 30,000 models for every generic model of human atrial electrophysiology. This candidate population is then simulated in conditions closely resembling experiment, and a subset of models with biomarkers within experimentally observed ranges are retained for further analysis. The SMC algorithm searches parameter space until it finds a certain number of parameter sets (here 800) that yield in silico models whose outputs agree with our ex vivo AP biomarkers. To locate those target parameter sets, SMC repeatedly resamples models whose AP biomarkers are furthest away from the desired values and then applies Markov chain Monte Carlo steps to these models to explore the parameter space and ensure that most models in the population remain unique. The jumping distribution used for these Markov chain Monte Carlo steps is a three-component Gaussian mixture fitted to the locations of all models in the parameter space. A detailed description of the SMC method is provided in SMC Sampling for Experimentally Calibrated Populations of Models in the APPENDIX.
Calibration criteria from ex vivo recordings. Calibration criteria involved retaining the models whose AP and/or ionic properties were in range with ex vivo recordings. Most studies with the use of the experimentally calibrated population-of-models methodology had access to experimental data on AP biomarkers only (25). To investigate the importance of information contained within AP biomarkers versus current densities, we first constrained candidate populations using solely AP data. The experimentally derived constraints imposed on in silico APs consisted of the following: 1) minimum/maximum values of APD20, APD50, APD90, RMP, APA, and the difference APD90 Ϫ APD50 (triangulation) at all pacing frequencies ( Table 2); 2) constraint of linearity between APD20 and APD50 observed when ex vivo biomarkers were investigated for covariation (Fig. 2); 3) APD90 rate dependence: the values of APD90 at 3 and 0.25 Hz generated by each in silico model were required to fulfill the inequality 1.  A: AP biomarkers used in this study. APD20, APD50, and APD90, AP duration at 20%, 50%, and 90% repolarization, respectively (in ms); RMP, resting membrane potential (in mV); APA, action potential amplitude (in mV); V, voltage. B-E: outputs of generic Courtemanche, Maleckar, and Grandi models compared with experimental AP traces (B), AP biomarkers (C), and ionic current densities (D and E). In C, outputs of Courtemanche, Maleckar, and Grandi models are represented with a diamond, square, and circle, with ex vivo data plotted as crosses. E: two peaks in L-type Ca 2ϩ current (ICaL) density were observed ex vivo. To capture these, modifications to half-potential of activation (Vact) of the ICaL d-gate were introduced to generic models (see Modifications to the generic models to capture two peaks in ICaL density in the APPENDIX). whose APD90 differed by Ͼ5 ms between two final beats were rejected, as in Ref. 52. In the case of LHS, all experimental constraints on the AP properties were applied at once. Under SMC, the calibration criteria in points 1 and 2 were applied sequentially for every biomarker (SMC Sampling for Experimentally Calibrated Populations of Models in the APPENDIX); to reduce computational burden, criteria 3 and 4 were applied after the SMC algorithm stopped, i.e., after generating 800 models that fulfilled constraints in points 1 and 2.
In silico stimulation protocols. In silico populations were simulated to elicit cellular APs and I to, IKur, ICaL, and IK1 densities in conditions resembling experiments as closely as possible (for details, see Computational stimulation protocols and Numerical methods and data analysis in the APPENDIX). Briefly, we assumed that ionic concentrations in the pipette and extracellular solutions used in experiments corresponded to the intracellular and extracellular concentrations in the models. The intracellular Na ϩ and K ϩ and extracellular Na ϩ , K ϩ , and Ca 2ϩ concentrations were held constant in time and to their experimental equivalents. The temperature was set to 310.15 K in the simulations where the APs and Ito, IKur, and ICaL current densities were elicited, whereas in the simulations of IK1 density, the temperature was 295.15 K. To elicit APs, each model was stimulated for 100 beats at every frequency, following the framework designed in Ref. 52. The duration and amplitude of stimulus current were matched to the median experimental values. In voltage-clamp simulations, pharmacological inhibition of specific ion channels was modeled by setting the corresponding maximal conductances to zero.

Populations of models with wide variations in ionic densities
can yield APs overlapping with experiments, even if generic models are initially far away from the specific experimental cohort. Figure 3 shows the AP biomarkers and ionic densities in Maleckar-based and Grandi-based in silico populations produced with the LHS algorithm following calibration with ex vivo AP data. In both populations, substantial variability in maximal conductances was sufficient to generate a wide range of AP characteristics that overlapped with our ex vivo data set, even if the generic model outputs were far from the ex vivo data under consideration (Fig. 1). In the Maleckar-based population, calibration with ex vivo AP biomarkers yielded n ϭ 1,493 in silico models with AP properties spanning full ranges of experimentally observed variability at the AP level (Fig. 3,  A and B). Simulated ionic densities largely encompassed ex vivo recordings of inward Ca 2ϩ and outward K ϩ currents, with in silico I K1 and I CaL spanning a wider range of values than in our experiments (Fig. 3C).
A similar picture could be seen in the Grandi-based population. Of the N ϭ 30,000 candidate models generated by LHS, n ϭ 554 produced APD and APA properties in range with experimental observations (Fig. 3, D and E), with simulated RMP biomarkers occupying the bottom half of the experimentally observed range of values. Ionic current densities in the Grandi-based population overlapped with ex vivo recordings, with some models exceeding peak I CaL density recorded in experiments (Fig. 3F).
Not all models succeeded in reproducing the variability of the AP recordings. Despite the wide variability in maximal conductances, the Courtemanche-based in silico population of n ϭ 1,553 models failed to cover experimental AP ranges, particularly the early repolarization phase (Fig. 4). This is because of the spike-and-dome morphology of the population's APs, capable of producing APD 20 and APD 50 only at the low Table 2. Values are ranges with means Ϯ SD in parentheses. APD20, APD50, and APD90, action potential duration at 20%, 50%, and 90% repolarization, respectively; RMP, resting membrane potential; APA, action potential amplitude; Tri90-50, triangulation.

H898 FROM IONIC TO CELLULAR VARIABILITY IN HUMAN ATRIAL MYOCYTES
end of values recorded in our ex vivo data set ( Fig. 4), thus failing to yield the long APD 20 and APD 50 seen in our experimental data set. The findings of the Courtemanchebased, Maleckar-based, and Grandi-based populations were replicated with the SMC algorithm (Figs. 5-7), highlighting the independence of our results from the specific method used to randomly sample the conductance values in the population. In the remainder of this report, we focus our attention on the Maleckar-based and Grandi-based in silico populations, since they captured experimental ranges of variability in the AP.
Calibration with AP data constrains currents affecting upstroke and resting potential, but current redundancy in repolarization allows wide ranges of variability in currents impacting APD. Figure 8 shows scatterplots of key maximal conductances, and their correlation with AP biomarkers for the Maleckar-based and Grandi-based populations, calibrated with AP data. In the Maleckar-based population, models with Na ϩ conductance below Ϫ59% of baseline are absent from the population because of the critical role of I Na in the generation of the AP amplitude (Fig. 8B). Models with low Na ϩ -K ϩ pump conductance (G NaK ) and low inward rectifier conductance (G K1 ) are absent from both AP-calibrated populations because of the critical importance of I NaK and I K1 for RMP (Fig. 8, B and D). The current conductances regulating cellular repolar-ization are unconstrained by calibration with AP data because multiple ionic currents regulate APD biomarkers (for instance, APD 20 is determined by a balance of the following four ionic currents: I to , I Kur , I CaL , and I Na ; Fig. 8

, B and D).
We can now rationalize the impact of AP-level calibration on the simulated current densities shown in Fig. 3, C and F. It is important to note that the experimental measurements of a specific current density and cellular AP were not obtained simultaneously from the same cell, i.e., our experimental recordings of I to , I Kur , I K1 , and I CaL do not underpin our experimental APs (furthermore, the cells used for measurements of one current were different from those used for recordings of another current). Figure 3, C and F, highlights large variability in peak current density observed in experiments and simulations. In both Maleckar-based and Grandi-based populations calibrated with AP-level data, the ranges of ultrarapid K ϩ conductance (G Kur ) and transient outward K ϩ conductance parameters spanned Ϯ100% of their baseline values (Fig. 8, A and C), translating into simulated I Kur and I to densities that cover the bulk of our experimental data (Fig. 3, C and F).
Differences are observed in the ranges of I K1 density between Maleckar-based and Grandi-based populations. In both populations, high G NaK can compensate for low G K1 and vice versa, as shown in Fig. 8, A and C. In the Maleckar-based Linear relationship between action potential (AP) duration (APD) at 20% and 50% repolarization (APD20 and APD50) properties at five pacing frequencies ex vivo. Experimental data are shown as crosses, dashed-dotted lines are the least-squares fits to the data (corresponding R 2 value given at the top of each plot), and solid black lines represent constraints corresponding to aAPD20 ϩ bupper Ն APD50 Ն aAPD20 ϩ blower imposed on APD20 and APD50 properties generated by in silico models. The constant a was obtained from the least-squares fit, whereas blower and bupper were determined from the values of the outermost experimental data points. population, some models displayed a peak I K1 density of up to Ϫ24.62 pA/pF (Fig. 3C), exceeding the value of Ϫ10.33 pA/pF seen in our ex vivo data, whereas the Grandi-based population produced I K1 density in range with our experimental recordings (Fig. 3F). This interpopulation difference is the result of distinct experimental data sets used in the construction of the generic Maleckar and Grandi models. The generic Maleckar model yielded I K1 density of Ϫ12.31 pA/pF (at Ϫ100 mV) and was based on previous experimental data (3) where I K1 density of Ϫ9.83 Ϯ 0.96 pA/pF (at Ϫ90 mV) was reported. The Grandi model, derived from a model of a human ventricular myocyte, takes into account the sixfold smaller density of I K1 in atria relative to ventricles (15,42) and yields a baseline I K1 density of Ϫ5.49 pA/pF, close to the experimentally measured value of Ϫ3.5 Ϯ 0.3 (at Ϫ120 mV) previously reported (55). The threefold difference in experimental measurements of I K1 density previously reported (3,55) provides another illustration of the wide heterogeneity in the human population.
In both in silico populations calibrated with AP data, I CaL was not the main determinant of AP, as evidenced by its weak correlation with AP properties (Fig. 8, B and D). Accordingly, the L-type Ca 2ϩ conductance (G CaL ) parameter spanned Ϯ100% of generic model values in both populations (Fig. 8, A and C). In some of the Maleckar-based and Grandi-based models, I CaL density was as high as Ϫ7 and Ϫ10 pA/pF, both of which exceeded the maximal I CaL density of Ϫ4.35 pA/pF in our ex vivo recordings (Fig. 3, C and F). However, peak I CaL density ranging from Ϫ0.1 to Ϫ9 pA/pF (see Fig. 1 in Ref. 8) and from~0 to Ϫ55 pA/pF (see Fig. 3 in Ref. 49) have been reported in patient cohorts in sinus rhythm similar to ours. Therefore, models with larger peak I CaL are still representative of the variability observed in sinus rhythm.
Further calibration with ionic currents active during repolarization constrains current density but not cellular phenotypic variability in AP. We proceeded to further calibrate the two in silico populations using voltage-clamp I CaL and I K1 measurements by retaining models whose peak current densities do not exceed the maximal values in our ex vivo data set. Figures 9 and 10 show the impact of this additional calibration , ultrarapid K ϩ current (IKur), inward rectifier (IK1), and L-type Ca 2ϩ current (ICaL) (with peaks at 10 and 20 mV (C) in the Maleckar-based population of n ϭ 1,493 models following calibration with AP biomarker data compared with ex vivo measurements. D-F: analogous plots for the Grandi-based in silico population of n ϭ 554 models. Both in silico populations were generated with Latin hypercube sampling. Throughout, in silico "accepted" models, meaning those models that fulfill experimental AP calibration constraints, are plotted in dark gray; A and D additionally show AP biomarkers of in silico "rejected" models in light gray. Ex vivo data are shown as crosses. RMP, resting membrane potential; APA, AP amplitude; APD20, APD50, and APD90, AP duration at 20%, 50%, and 90% repolarization. on maximal conductances, AP biomarkers, and peak current densities in the Maleckar-based and Grandi-based populations, respectively. In both populations, additional calibration constraint with I CaL density had no effect on maximal conductances other than G CaL or on AP biomarkers (Figs. 9A, 10A, 11, and 12). This is in line with the weak correlation of G CaL with AP properties in Fig. 8 and further illustrates the limited role of I CaL on cellular phenotypic variability in the AP in our data set.
In contrast, the impact of additional calibration with peak I K1 density differed between in silico populations. The Maleckar-based population calibrated with AP data contained many models with I K1 current exceeding the peak value seen in our ex vivo recordings. Thus, additional calibration on I K1 density substantially reduced the number of viable models and impacted not only the distribution of G K1 but also of G NaK (since those parameters are correlated; Fig. 8A). Because of the I K1 importance in determining APD 90 , few models with short APD 90 survived this additional calibration step (Fig. 9B). In contrast, in the Grandibased population calibrated with AP-level data, simulated I K1 density effectively spanned the ex vivo range; hence, the calibration step on peak I K1 density had little impact on the parameter space or AP biomarkers (Fig. 10). Following the calibration with both AP and voltage-clamp level information, we obtained Maleckar-based and Grandi-based in silico populations with AP and ionic densities in range with our ex vivo data set (Fig. 9, C and D, and Fig. 10, C and D). following calibration with AP biomarker data compared with ex vivo measurements. In silico population was generated with Latin hypercube sampling. In silico accepted models, i.e., those models that fulfill experimental AP calibration constraints, are plotted in dark gray; A additionally shows AP biomarkers of in silico rejected models in light gray. Ex vivo data are shown as crosses. , and ionic current densities (C) of transient outward K ϩ current (Ito), atrial-specific ultrarapid K ϩ current (IKur), inward rectifier K ϩ current (IK1), and L-type Ca 2ϩ current (ICaL) (with peaks at 10 and 20 mV) in the Maleckar-based population of n ϭ 754 models following calibration with AP biomarker data compared with ex vivo measurements. The in silico population was generated with sequential Monte Carlo. In silico accepted models, i.e., those models that fulfill experimental AP calibration constraints, are plotted in gray dots, whereas ex vivo data are shown as crosses.  , and ionic current densities (C) of transient outward K ϩ current (Ito), atrial-specific ultrarapid K ϩ current (IKur), inward rectifier K ϩ current (IK1), and L-type Ca 2ϩ current (ICaL) (with peaks at 10 and 20 mV) in the Grandi-based population of n ϭ 479 models following calibration with AP biomarker data compared with ex vivo measurements. In silico population was generated with sequential Monte Carlo. In silico accepted models, i.e., those models that fulfill experimental AP calibration constraints, are plotted in light gray, whereas ex vivo data are shown as crosses. To further explore variability in human atrial models, we analyzed Ca 2ϩ transient (CaT) properties in both in silico populations. Figure 13 shows the main CaT types in both populations calibrated with experimental AP and voltageclamp information, whereas Fig. 14 shows the impact of calibration with peak I CaL density and highlights the ionic determinants of CaTs. Given that we considered Ϯ100% variation in many ionic densities, including in the sarcoplasmic reticulum release and uptake currents (J rel and J up ; Figs. 11 and 12), a wide variability in the biomarkers and morphology of the simulated CaTs can be expected.
The models within the Maleckar-based population exhibited the following three CaT types (Fig. 13, A and B): the morpho-logically normal transient, seen in 60% of the models; Ca 2ϩ oscillations in diastole, observed in 21% of the models and previously reported in ex vivo multicellular tissue preparations (33); and spontaneous Ca 2ϩ release from the sarcoplasmic reticulum during diastole, observed in 19% of the models. The proportions of models with the three CaT types did not change substantially with I CaL calibration (Fig. 14A), illustrating that the CaT phenotypes in the Maleckar-based population are independent of I CaL density.
The Grandi-based population calibrated with AP and voltage-clamp level data contained three main CaT types (Fig. 13, D and E) as follows: the morphologically normal transient, exhibited by 55% of the models; the spike-and-dome CaT characterized by the presence of two peaks in CaT curve during systole and seen in 42% of the models; and spontaneous Ca 2ϩ Fig. 7. In silico action potential (AP) biomarkers at 0.25, 1, and 3 Hz (A) and full AP traces (B) in the Courtemanche-based population of n ϭ 819 models following calibration with AP biomarker data compared with ex vivo measurements. In silico population was generated with sequential Monte Carlo. In silico accepted models, i.e., those models that fulfill experimental AP calibration constraints, are plotted as light gray, whereas ex vivo data are shown as crosses. release in diastole, observed in 2% of the models. The proportion of models with the spike-and-dome CaT morphology increased significantly with I CaL calibration (Fig. 14C), suggesting that low I CaL may promote this CaT phenotype in Grandi-based models.
In both in silico populations, fluctuations in diastolic Ca 2ϩ (inclusive of Ca 2ϩ oscillations in diastole and spontaneous Ca 2ϩ release phenotypes) were observed in models with low NCX conductance (G NCX ; Fig. 13, C and F). Additionally, plots of ionic currents involved in CaT generation implicated a release from the sarcoplasmic reticulum during diastole (Figs.  15 and 16). In a study on guinea pig hearts ex vivo, Plummer et al. (33) reported multicellular diastolic Ca 2ϩ fluctuations, observing that RyR inhibition abolished those fluctuations, whereas increasing RyR open probability exacerbated them. Consistent with this, a reduction of J rel by 90% suppressed most diastolic Ca 2ϩ fluctuations, whereas increasing J rel by 100% had the opposite effect in both in silico populations (Tables 3 and 4). Separately, increasing G NCX to normal levels abolished most Ca 2ϩ fluctuations in diastole in both populations (Tables 3 and 4).
An interesting feature of some models within the Grandibased population is the spike-and-dome CaT morphology, coinciding with strong extrusion of Ca 2ϩ from the myocyte (high G NCX ; Fig. 13F) and a weak injection of Ca 2ϩ in the cytosol [low J rel (Fig. 13F) and/or low I CaL (Fig. 14C)]. The spike in the CaT during systole coincides with a prominent "notch" in forward-mode I NCX (Fig. 16). Accordingly, 50% G NCX block to diminish strong extrusion of Ca 2ϩ or a 50% J rel /I CaL increase for a stronger injection of Ca 2ϩ in the cytosol (to compensate for the effect of I NCX notch), all substantially reduced the number of models with this CaT phenotype (Table  5). In summary, in the Maleckar-based population, the biomarkers of morphologically normal CaTs are mainly deter- mined by sarcoplasmic reticulum release and uptake currents (J rel and J up ; Fig. 14B), whereas G NCX and J rel control the appearance of diastolic Ca 2ϩ fluctuations. In the Grandi-based population, G NCX and J up are responsible for the properties of morphologically normal CaTs (Fig. 14D); additionally, G NCX and J rel control both the diastolic Ca 2ϩ fluctuations and spikeand-dome CaT phenotypes, with the latter also modulated by G CaL .

DISCUSSION
In this report, we describe a combined experimental and computational investigation to deepen our quantitative under-standing of cellular and ionic variability in human right atrial myocytes. We developed experimentally calibrated populations of in silico models based on an ex vivo data set of AP and current density recordings obtained in right atrial myocytes from patients undergoing coronary revascularization or aortic valve replacement. Our main findings are as follows: 1) variability in human atrial APs and ionic currents is wide, with large differences observed a) within and between experimental data sets from different patient cohorts (8,39,49,50) and b) between APs/ionic densities in experimental recordings versus generic in silico models; 2) ex vivo cell-to-cell variability in peak current densities from human atrial cardiomyocytes is substantial and can translate to a Ϯ100% variation in maximal conductances in populations of in silico models of human electrophysiology; 3) populations of models with a substantial variability in maximal conductances can yield AP characteristics overlapping with ex vivo data set even if outputs of in silico generic models are far away from experiment; 4) calibration with AP data constrains currents affecting upstroke and resting potential, but current redundancy in repolarization allows wide ranges of variability in currents impacting APD; 5) additional calibration with ionic currents active during repolarization constrains current density but not cellular phenotypic variability in AP; and 6) experimentally calibrated populations of atrial electrophysiological models exhibit three Ca 2ϩ transient phenotypes, mainly determined by NCX and RyR/ SERCA magnitudes rather than I CaL . The appearance of Ca 2ϩ fluctuations in diastole can be caused by low I NCX , in addition to the previously reported effect of elevated sarcoplasmic release (33).
Our experimental data yielded relatively short APD biomarkers (e.g., APD 90 ranged between 69 and 143 ms at 1 Hz; Table 2) and APs with a triangular morphology. These data are consistent with findings reported by other groups in atrial myocytes from patients in sinus rhythm (41,45,55,56). Other studies have reported APD 90 values between 190 and 440 ms at 1 Hz in atrial trabecular preparations from patients in sinus rhythm (39), with both triangular and spike-and-dome AP shapes observed. Methodological aspects (e.g., solution composition, temperature, cellular vs. multicellular recordings) and differences in patient cohorts (e.g., mitral valve vs. coronary revasculization patients) are likely to underlie these discrepancies. Such differences also extend to generic in silico models, depending on the data on which they were based. Given the specific electrophysiological characteristics of our cohort, our results might therefore be only generalizable in the case of short triangular APs in sinus rhythm, as also reported by others (41,45,55,56).
At the level of current densities, we demonstrated that variability in cellular ionic profiles ex vivo is substantial and can translate to a Ϯ100% variation in maximal conductances in the populations of in silico models. The wide range of variability supports the hypothesis that biological systems are generally "sloppy" (17), i.e., tolerant to significant variations of many parameter combinations. We investigated the ability of Courtemanche, Maleckar, and Grandi human single cell models to capture the variability observed in our ex vivo data set. All three generic models yielded in silico populations that spanned the ranges of ex vivo variability in APD 90 . However, because of the spike-and-dome morphology of its APs, the Courtemanche population was limited in producing models with long APD 20 and APD 50 as seen in the experimental data set and captured successfully by Maleckar-and Grandi-based populations. Despite the wide range of variability in maximal conductances, the resting potential spanned the bottom onehalf of experimentally measured values in all in silico populations. Aside from I K1 and I NaK , RMP is critically determined by the concentrations of ionic species on either side of the cell membrane (predominantly K ϩ ). Introducing a likely experimental variability in K ϩ concentrations could produce models with higher RMP, as previously shown (26,50). Additional sources of variability could also include cell-to-cell differences in cell volume and membrane capacitance values, as previously shown (34).
The experimental data set used in this work was particularly restrictive in terms of the short APs observed in the investigated cohort. This was in fact reflected when the in silico populations were calibrated against AP biomarkers, with a rate of model pruning significantly higher than in previous studies (4,20). Our present results further demonstrate that calibration against AP biomarkers can tightly restrict variability in ion channel densities when these strongly correlate with specific AP biomarkers (such as upstroke and resting potential). Contrastingly, redundancy in plateau and repolarization currents allows wider ranges of variability, and additional ionic density measurements may prove useful in further constraining some of these currents (as with I CaL here in both Maleckar-and Grandi-based populations). Although of relevance for previous studies on cell-to-cell variability in which only AP recordings were available and similar ranges of variation in maximal conductances were explored, it is important to recall that those studies aimed to investigate larger, more heterogeneous patient cohorts (20,39). Given that the mean current density of plateau and repolarization currents has been shown to vary more than fourfold in sinus rhythm in humans (8,49), larger I CaL than in our data may exist at the larger population scale.
A limitation of our experimental data set is that cellular APs and current densities were obtained in different cells. Nevertheless, cells were randomly used from all patients for both AP and ionic current measurements. This should prevent bias and provide a robust rationale for combining both types of measurements. It should be noted that obtaining all measurements Fig. 11. Box plots of maximal current conductances in the Maleckar-based in silico population of models following calibration with AP, AP ϩ ICaL, and AP ϩ ICaL ϩ IK1 data, shown in light gray, dark gray, and black, respectively. Black dots indicate outliers. Fig. 12. Box plots of maximal current conductances in the Grandi-based in silico population of models following calibration with action potential (AP), AP ϩ L-type Ca 2ϩ current (ICaL), and AP ϩ ICaL ϩ inward rectifier K ϩ current (IK1) data, shown in light gray, dark gray, and black, respectively. Black dots indicate outliers.
in the same cell is not practical, since cells are unlikely to survive exposure to both voltage-and current-clamp configurations, different electrode and bath solutions, and multiple pharmacological interventions and washout periods. It is, however, important to stress that current densities are dynamically changing quantities in vivo, and experimental investigations ex vivo offer valuable but limited snapshots. Considering these complexities, the goal of the experimentally calibrated populations of models is indeed to explore many ionic profiles yielding viable APs by generating larger pools of cellular models than the number of cells that is experimentally feasible to investigate, to overcome some of the experimental limitations.
Further investigations using the resulting in silico populations were conducted to analyze their underpinning Ca 2ϩ transients and, in particular, to investigate mechanisms of diastolic Ca 2ϩ fluctuations. These were mainly determined by NCX and RyR/SERCA magnitudes rather than I CaL . Fluctuations were more frequent in the Maleckar-based than Grandibased population, in some cases leading to delayed afterdepolarizations, consistent with previous findings (6). Importantly, fluctuations in diastolic Ca 2ϩ correlated with low NCX expression in both populations, indicating model independence. An implication is that NCX regulates the appearance of diastolic Ca 2ϩ disturbances in human atrial cardiomyocytes, in addition to previously reported effects of RyRs (33). In summary, our study provides further experimental and computational evidence on the large variability in ionic densities and AP recordings in human atrial myocytes and their potential implications in explaining cellular phenotypic variability in transmembrane voltage and Ca 2ϩ transient dynamics.
Integration between experiments and simulations extends the existing knowledge on the ionic and subcellular mechanisms underlying variability in the human atrial AP and their possible consequences in intracellular Ca 2ϩ dysregulation.

Experimental Data Set
Human atrial myocytes were isolated from the right atrial appendages of patients in sinus rhythm undergoing coronary artery bypass graft/aortic valve replacement/mitral valve replacement in Oxford, UK, as previously described (36). Table 1 shows detailed information on the patient cohort used to perform whole cell voltage and currentclamp experiments in this study.
AP measurements. Patch-clamp experiments were performed with the pipette solution containing (in mM) 130 potassium aspartate, 10 HEPES, 5 Na 2ATP, 2 MgCl2, 5 CaCl2, and 11 EGTA (pH adjusted to 7.2 using KOH), whereas the extracellular solution contained (in mM) 140 NaCl, 5.4 KCl, 1.8 CaCl 2, 1.2 MgCl2, 10 glucose, and 5 HEPES (pH adjusted to 7.4 using NaOH). A preexperimental protocol was run to determine the stimulus intensity: APs were recorded in response to 0.25-to 1.5-nA depolarizing current pulses lasting 4 ms with an interval of 1 s using an Axopatch 200 B amplifier in current-clamp mode. The stimulus intensity could be 1.5 times the diastolic threshold. APs from a total of n ϭ 29 cells from N ϭ 8 patients were measured at five stimulation frequencies (0.25, 0.5, 1, 2, and 3 Hz). Every cell was stimulated a minimum of 25 times; the recorded APs were subsequently averaged to yield a single mean AP trace used to determine the AP properties for each cell. For each cellular AP, RMP, APA, APD 20, APD50, and APD90 were quantified.
Voltage-clamp data. Whole cell voltage-clamp experiments were used to determine I to, IKur, IK1, and ICaL, as previously detailed in Ref. 36.
To measure Ito, outward currents were evoked by steps in the range of Ϫ70 to ϩ50 mV for 700 ms from a holding potential of Ϫ80 mV (sampling rate: 5 kHz); a prestep of 200 ms to Ϫ40 mV was used to inactivate Na ϩ current. Ito was measured as the steady-state current subtracted from peak outward current at the start of the depolarizing step and normalized to membrane capacitance value. Steady-state currents were measured at the end of the depolarizing steps. I to was measured in n ϭ 30 cells from N ϭ 7 patients.
To measure I Kur, outward currents were evoked by steps in the range of Ϫ35 to ϩ45 mV for 500 ms from a holding potential of Ϫ80 mV; a prestep of 200 ms to Ϫ35 mV was used to inactivate Na ϩ current. Sustained current (Isus) was measured as the amplitude of the current remaining at the end of the pulse relative to the zero current level at the end of the prestep to Ϫ35 mV. 4-Aminopyridine (4-AP; 50 mol/l)-sensitive I Kur was obtained by subtracting Isus from the remaining current after 4-AP inhibition (53). IKur was measured in n ϭ 30 cells from N ϭ 7 patients.
In voltage-clamp experiments measuring I Kur and Ito, the pipette solution contained (in mM) 120 potassium aspartate, 20 KCl, 5 Na 2ATP, 2 MgCl2, 5 EGTA, and 10 HEPES (pH adjusted to 7.2 using KOH; pCa 7.0). The extracellular solution contained (in mM) 140 NaCl, 5.4 KCl, 1.8 CaCl 2, 1.2 MgCl2, 10 glucose, 5 HEPES, and 0.5 CdCl 2 (pH adjusted to 7.4 using NaOH). Cd 2ϩ was used in experiments to block L-type Ca 2ϩ channels. However, Cd 2ϩ are known to shift the steady-state activation and inactivation of K ϩ repolarization currents to more positive potentials by ϳ20 mV (46). We account for this by shifting experimental current-voltage curves of K ϩ currents by Ϫ20 mV. I K1 was obtained by digital subtraction of current traces, recorded with a ramp pulse from Ϫ100 to ϩ40 mV, in the presence or absence of 1 mmol/l BaCl 2 (53), as previously described (51). The extracel-  ICaL was elicited with step depolarization protocols using test potentials in the range of Ϫ50 to ϩ50 mV in 10-mV increments with a holding potential at Ϫ40 mV (43). Recordings were begun and finished after the Ca 2ϩ current amplitudes and access resistance had stabilized. ICaL was measured as the difference between the current at the end of 250-ms test pulse and the peak inward current. ICaL was measured in n ϭ 18 cells from N ϭ 6 patients. The cells could be split into two groups based on the value of the transmembrane voltage at which the peak inward current was detected, since peak I CaL was registered at 10 and 20 mV in 8 and 10 cells, respectively. Pipette solution contained (in mM) 120 Cs-Aspartate, 10 TEA-Cl, 5 MgATP, 2 MgCl 2, 1 CaCl2, 11 EGTA, and 10 HEPES (pH 7.2 with CsOH). The extracellular solution contained (in mM) 140 TEA-Cl, 5.4 CsCl, 1.2 MgCl2, 5 HEPES, 11 glucose, and 1.8 CaCl2 (pH 7.4 with CsOH).

Generic Models of Human Atrial Electrophysiology
Three recent models of human atrial electrophysiology by Courtemanche et al. (7), Maleckar et al. (23), and Grandi et al. (15) were used as a basis for the experimentally calibrated atrial cell populations. Throughout the report, we refer to these models by their first authors' names for simplicity.
In brief, the three models include the major transmembrane currents responsible for AP generation such as I Na, IKur, Ito, IKr, IKs, ICaL, IK1, INaK, and INCX. The models also incorporate a description of intracellular Ca 2ϩ handling that includes the uptake and release currents Jup (SERCA) and Jrel (RyR) of the sarcoplasmic reticulum as well as K ϩ , Na ϩ and Ca 2ϩ homeostasis regulating intracellular ionic concentrations.
However, the three models differ in the mathematical formulations of most of the transmembrane currents, the intracellular Ca 2ϩ handling subsystem, and the number of model compartments. The intracellular Ca 2ϩ handling subsystem in the Courtemanche model is based on the guinea pig ventricular model of (22), whereas the Grandi model uses the formulation of the rabbit ventricular Ca 2ϩ handling previously described (44). In the Maleckar model, the Ca 2ϩ subsystem is described as in the rabbit atrial model previously published (21). Intermodel differences related to the number of cellular compartments are also seen. The Courtemanche model is simplest in this regard, consisting of the extracellular and intracellular spaces between which the transmembrane currents flow. The Maleckar model additionally comprises the cleft space outside of the sarcolemma, where local accumulation and depletion of ions is permitted. This cleft space acts as a buffer between the intracellular and extracellular ionic solutions and can be thought of as an effective extracellular space. The Grandi model consists of the intracellular, subsarcolemmal, and extracellular spaces, with the subsarcolemmal space bridging the intracellular and extracellular compartments and acting as an effective intracellular space in terms of ionic concentrations seen by transmembrane currents.
As a result of these differences, the Courtemanche model generates a spike-and-dome AP morphology with a slowly decaying intracellular CaT, whereas the Maleckar model produces a triangular AP with a larger CaT amplitude and faster decay. The Grandi model produces an AP that is triangular but lower in amplitude than in the other models, and its CaT is slower to rise and decay.
Modifications to the generic models to capture two peaks in ICaL density. Figure 1, D and E, shows a comparison of Ito, IKur, IK1, and I CaL densities recorded ex vivo with those yielded by the three generic in silico models. Simulated current-voltage curves of K ϩ repolarization currents in all three models qualitatively follow experimental recordings. In silico, ICaL density is maximized at 10 mV in the Courtemanche and Maleckar models and at 0 mV in the generic Grandi model (7,15,23). Analysis of ex vivo ICaL voltage-clamp data revealed two groups of cells with I CaL density maximized at 10 and 20 mV, respectively (Fig. 1E). To capture experimental variability in I CaL peak density location, L-type Ca 2ϩ channel gating was modified in all three generic models. Parameters that control the membrane potential maximizing ICaL density include half-potentials of activation (Vact) and inactivation (V inact). In the generic Courtemanche, Maleckar, and

SMC Sampling for Experimentally Calibrated Populations of Models
The SMC approach imposes the experimentally motivated calibration constraints in distinct stages (this contrasts with LHS, where all experimentally motivated constraints can be applied simultaneously to all models within the candidate population). At each stage, SMC requires that the values of a specific AP biomarker produced by all models in the population fall within the experimentally permitted ranges for all pacing frequencies. Introduction of a new set of constraints can potentially result in sharp reductions in the regions of the parameter space that are viable. To "soften" this effect, constraints are not introduced in a strict fashion.
Expressed in simplified form, the algorithm proceeds as follows: 1) initializes particles randomly throughout the search region of the parameter space (base values of model parameters Ϯ 100%); 2) tests each particle for how close it falls to the acceptable range of values for the first biomarker at all pacing frequencies (this distance is the 2-norm of the vector of distances away from the acceptable range for each biomarker at each pacing frequency); 3) takes all particles that fully satisfy the constraint for all pacing frequencies (distance of 0) and then additional particles selected per increasing distance until at least one-half of the particles have been chosen (the distance from the acceptable range of values of the furthest particle becomes the new acceptable limit); 4) copies the remaining particles (those outside the new acceptable limit) on random members of the set of particles that were chosen; 5) performs several iterations (here 5) of Markov chain Monte Carlo on the copied particles to move them to unique locations in the parameter space that are still within the currently acceptable distance away from the constraint; 6) repeats steps 2 to 5 until all particles correspond to models that produce values for the current biomarker that fall within the desired ranges for each pacing frequency; and 7) repeats steps 2 to 6 for each biomarker, ensuring that, during trialed updates, no particles move to locations that would break any of the previously satisfied constraints.
The method is explained in full mathematical detail (9). Following the approach described therein, the Markov chain Monte Carlo steps are performed using a jumping distribution constructed by fitting a three-component Gaussian mixture model to the set of particles selected before each copy step.

Computational Stimulation Protocols
Computational AP protocol. We assumed that ionic concentrations in the pipette and extracellular solutions used in experiments corresponded to the intracellular and extracellular concentrations in the models. The intracellular Na ϩ and K ϩ and extracellular Na ϩ , K ϩ , and Ca 2ϩ concentrations were held constant in time. Each model was stimulated for 100 beats at every frequency. Stimulus current duration and amplitude were adjusted in the models to capture the typical amount of charge injected in the cells ex vivo (60 pC). In the simulations with the candidate population based on the Maleckar model, the intracellular K ϩ , Na ϩ , and Mg 2ϩ and effective extracellular K ϩ , Na ϩ , and Ca 2ϩ concentrations were set to 130, 10, 2, 5.4, 140, and 1.8 mM, respectively. In the simulations with the populations based on the Grandi model, the effective intracellular K ϩ , Na ϩ , Cl Ϫ , and Mg 2ϩ and extracellular K ϩ , Na ϩ , Ca 2ϩ , and Cl Ϫ concentrations were set to 130, 10, 14, 2, 5.4, 140, 1.8, and 151.4 mM, respectively. For the candidate populations derived from the Courtemanche model, the intracellular K ϩ and Na ϩ and extracellular K ϩ , Na ϩ , and Ca 2ϩ concentrations were set to 130, 10, 5.4, 140, and 1.8 mM, respectively. In an attempt to better mimic experimental conditions ex vivo, the (effective) intracellular Na ϩ and K ϩ and (effective) extracellular Na ϩ , K ϩ , and Ca 2ϩ concentrations were held constant in time. In the candidate populations based on the Courtemanche model, this meant clamping the intracellular K ϩ and Na ϩ (the model lacks the effective extracellular space of the Maleckar model or the effective intracellular space of the Grandi model). In the populations based on the Maleckar model, the effective extracellular concentrations of K ϩ , Na ϩ , and Ca 2ϩ and the intracellular concentrations of K ϩ and Na ϩ were clamped. In the candidate populations based on the Grandi model, the Na ϩ concentrations in the intracellular, effective intracellular, and cleft spaces (the latter being a part of the Ca 2ϩ -handling subsystem) were clamped.
Computational voltage-clamp protocols. In the voltage-clamp simulation protocol designed to elucidate Ito, the holding potential of Ϫ80 mV was maintained for 200 ms before the prestep to Ϫ40 mV for another 200 ms. After this, the transmembrane potential was set to the desired value (between Ϫ70 mV and 50 mV in 10-mV increments, with the exception of the increment between Ϫ10 and 10 mV, which was set to 0.0001 mV) for 700 ms followed by the voltage step to Ϫ40 mV for 200 ms and then to Ϫ80 mV for another 200 ms.
In the voltage-clamp simulation protocol designed to elucidate IKur, the holding potential was set to Ϫ80 mV for 200 ms followed by a prestep to Ϫ35 mV for another 200 ms and then by the desired voltage step (between Ϫ35 and 45 mV in 10-mV increments) for 500 ms. After this, the membrane potential was set to Ϫ35 mV and then to Ϫ80 mV for 200 ms each. As in experiments, Isus was measured as the amplitude of the current remaining at the end of the pulse relative to the zero current level at the end of the prestep to Ϫ35 mV. To mimic the block of IKur by 4-AP, values of GKur were set to zero in the relevant simulations, and IKur was determined as the difference in the total transmembrane current Isus elucidated in the absence or presence of the simulated IKur block by 4-AP.
In the simulations of Ito and IKur, the inhibitory effect of Cd 2ϩ on ICaL was modeled by setting GCaL ϭ 0 (32). The (effective) intracellular K ϩ , Na ϩ , Mg 2ϩ , and Cl Ϫ concentrations were set to 140, 10, 2, and 24 mM, whereas the (effective) extracellular K ϩ , Na ϩ , Ca 2ϩ , and AP, action potential; ICaL, L-type Ca 2ϩ current; IK1, inward rectifier K ϩ current; Jrel, release current; INCX, Na ϩ /Ca 2ϩ exchanger current; NA, not applicable. 50% block of INCX abolishes the spike-and-dome Ca 2ϩ transient morphology in 55% of the models, whereas 50% increase in Jrel and ICaL abolish the spike-and-dome morphology in 30% and 51% of the models.
Cl Ϫ concentrations were set to 5.4, 140, 1.8, and 152.4 mM, respectively, where applicable. The temperature was set to 310.15 K.
In the protocols for measuring ICaL, the holding potential was set at Ϫ40 mV for 125 ms followed by a voltage step to the desired voltage (between Ϫ50 and 50 mV in 10-mV increments, with the exception of the increment between Ϫ10 and 10 mV, which was set to 0.0001 mV) for 250 ms. After this, the transmembrane potential was held at Ϫ40 mV for another 125 ms. The temperature was set to 310.15 K.
In the simulations of I CaL, all transmembrane K ϩ currents (including the background K ϩ current in the Courtemanche model and the junctional component of IKs in the Grandi model) were switched off by setting their conductances to zero to mimic experimental inhibition of K ϩ channels by cesium ions (32); Na ϩ currents (including the background currents) were also switched off because of the absence of Na ϩ in the pipette and extracellular solutions in ICaL voltage-clamp experiments. The absence of Na ϩ and K ϩ in experimental solutions was mimicked by setting the (effective) intracellular and extracellular Na ϩ and K ϩ concentrations to 0.0001 mM. The (effective) extracellular Ca 2ϩ and Cl Ϫ concentrations were set to 1.8 and 151.4 mM, whereas the (effective) intracellular Mg 2ϩ and Cl Ϫ concentrations were set to 7 and 16 mM. Exceptions here were the Grandi-based models, where intracellular Mg 2ϩ was set to 5.5 mM, since higher values of this concentration resulted in lack of numerical stability.
In the voltage-clamp simulation protocol designed to elucidate IK1, the holding potential was set to Ϫ80 mV for 50 ms followed by a prestep to Ϫ100 mV for another 50 ms and then by a voltage ramp from Ϫ100 to ϩ40 mV over the next 1,200 ms in fixed steps of 0.1167 mV/ms. After this, the membrane potential was set to Ϫ50 mV for 100 ms and then to Ϫ80 mV for 550 ms. The inhibitory effect of Ba 2ϩ on IK1 was mimicked by setting GK1 ϭ 0 in the relevant simulations, and IK1 was determined by subtraction of the total transmembrane current traces in the presence or absence of this simulated block. Temperature was set to 295.15 K. The (effective) extracellular concentrations of K ϩ , Na ϩ , Ca 2ϩ , and Cl Ϫ were set to 20, 120, 1.8, and 145.6 mM where applicable. The (effective) intracellular K ϩ , Na ϩ , Mg 2ϩ , and Cl Ϫ concentrations were set to 120, 8, 5, and 52 mM, respectively.
In all voltage-clamp simulations, the (effective) intracellular K ϩ and Na ϩ and (effective) extracellular K ϩ , Na ϩ , and Ca 2ϩ concentrations were held constant in time. The transmembrane potential was also held constant in time (i.e., dV/dt ϭ 0) with the exception of changes corresponding to the desired voltage steps. The stimulus current amplitude was set to zero in all simulations.

Numerical Methods and Data Analysis
The numerical simulations were performed using Chaste (31), an open-source software framework for modeling in computational biology. The CellML implementation of the single-cell human atrial electrophysiological model by Courtemanche et al. (7), Maleckar et al. (23), and Grandi et al. (15) was used. The ordinary differential equations were integrated with the software package CVODE (https:// computation.llnl.gov/casc/sundials/). A backward differentiation formula method used to integrate the differential equations had the absolute and relative error tolerances set to 10 Ϫ9 and 10 Ϫ10 . Electrophysiological properties were output every 1 ms. Postprocessing of AP traces to yield AP properties was performed in Chaste. Data analysis was performed with Matlab. Graphs were plotted using Matlab.