Disrupted Calcium Release as a Mechanism for Atrial Alternans Associated with Human Atrial Fibrillation

Atrial fibrillation (AF) is the most common cardiac arrhythmia, but our knowledge of the arrhythmogenic substrate is incomplete. Alternans, the beat-to-beat alternation in the shape of cardiac electrical signals, typically occurs at fast heart rates and leads to arrhythmia. However, atrial alternans have been observed at slower pacing rates in AF patients than in controls, suggesting that increased vulnerability to arrhythmia in AF patients may be due to the proarrythmic influence of alternans at these slower rates. As such, alternans may present a useful therapeutic target for the treatment and prevention of AF, but the mechanism underlying alternans occurrence in AF patients at heart rates near rest is unknown. The goal of this study was to determine how cellular changes that occur in human AF affect the appearance of alternans at heart rates near rest. To achieve this, we developed a computational model of human atrial tissue incorporating electrophysiological remodeling associated with chronic AF (cAF) and performed parameter sensitivity analysis of ionic model parameters to determine which cellular changes led to alternans. Of the 20 parameters tested, only decreasing the ryanodine receptor (RyR) inactivation rate constant (kiCa) produced action potential duration (APD) alternans seen clinically at slower pacing rates. Using single-cell clamps of voltage, fluxes, and state variables, we determined that alternans onset was Ca2+-driven rather than voltage-driven and occurred as a result of decreased RyR inactivation which led to increased steepness of the sarcoplasmic reticulum (SR) Ca2+ release slope. Iterated map analysis revealed that because SR Ca2+ uptake efficiency was much higher in control atrial cells than in cAF cells, drastic reductions in kiCa were required to produce alternans at comparable pacing rates in control atrial cells. These findings suggest that RyR kinetics may play a critical role in altered Ca2+ homeostasis which drives proarrhythmic APD alternans in patients with AF.


Intra-atrial heterogeneity and cAF remodeling
For the GPVm model, we implemented the changes described in Grandi et al. for the left atrium (LA) and right atrium (RA) [5]. In normal cells, the maximum conductance of IKur (gKur) is increased by 20% in the RA as compared to the LA. In cAF cells, the maximum conductances of IKur and Ito are differentially downregulated from their normal levels, with gKur decreased by 45% and 55% and gto decreased by 45% and 80% in the LA and RA, respectively. Additionally, modifications of the action potential under conditions of cAF included the following changes that occurred in both LA and RA: a 10% decrease in gNa, addition of a late component to the sodium current (INaL), a 50% decrease in gCaL, a 40% increase in IbarNCX, a 3-fold increase in the ryanodine receptor (RyR) Ca 2+ -dependent activation rate (koCa), a 25% increase in sarcoplasmic reticulum (SR) Ca 2+ leak, a 2-fold increase in gKs, and a 2-fold increase in gK1.

Sato-Bers RyR model implementation
Model parameters which were modified for implementation of the Sato-Bers RyR model [6] are listed in Table   2. The SR was divided into junctional (JSR) and network (NSR) compartments, the former of which contained the Ca 2+ buffer calsequestrin (CSQN). Equations for CSQN buffering in the Sato-Bers model were derived from Restrepo et al. [7]: The RyR equations were updated as described in Sato and Bers [6]: Equation S12 was modified to satisfy detailed balance.

Iterated map analysis
We used an iterated map analysis to derive Ca 2+ cycling stability criteria. For small SR load perturbations near steady state, total SR release ( ) and uptake ( ) on each beat changed linearly from beat to beat [8,9]: is the SR Ca 2+ release slope, is the SR Ca 2+ uptake factor, and and are the changes in total SR load and peak cytoplasmic Ca 2+ , respectively, from beats to . We did not consider Ca 2+ -induced Ca 2+ release (CICR) and SR leak separately since both were linearly dependent on SR load near steady state and their effects added linearly. Peak cytoplasmic Ca 2+ was defined as: (S15) where is the total Ca 2+ content in the cell at the start of beat . Equations S13-S15 were used to construct the first mapping equation describing the change in SR load [8][9][10]: (S16) Note that we did not assume to be zero. We also incorporated a new equation for net Ca 2+ efflux from the cell ( ), which depended linearly on peak cytoplasmic Ca 2+ for small perturbations near steady-state: (S17) is the sarcolemmal Ca 2+ efflux factor. Equations S13, S14, and S17 were used to find the linear least squares fit values of , , and based on , , , , and . 4 The Ca 2+ efflux term (Eq. S17) was used to construct the second mapping equation describing the change in total Ca 2+ content [10]: (S18) Though iterated map analysis lacking this second mapping equation has been previously used [9], the addition of the second equation provided more accurate theoretical predictions of Ca 2+ alternans thresholds in our simulations.
From Eq. S16 and S18, we obtained the following Jacobian matrix for the system [7,10]: (S19) The criteria for stability are that both eigenvalues of the Jacobian have absolute value less than 1: for Ca 2+ cycling to be stable.

Regression analysis
We used multivariable regression analysis methods from Sobie et al. to estimate the contribution of model parameters to the alternans threshold pacing cycle length (CL) [11]. Twenty model parameters (Table 1)  To determine the alternans threshold CL (output) for a given set of parameter scaling values (input), each cell was first paced to steady state at a CL of 400 ms. Then CL was progressively increased or decreased by 1 ms every 100 beats until APD alternans ceased (alternans ≤ 1%) or began (alternans > 1%), depending on whether alternans was present at a CL of 400 ms or not. Alternans threshold CL was defined as the shortest CL at which alternans did not occur. Any cell in which alternans persisted at CLs up to 750 ms or in which alternans was absent at CLs down to 100 ms was excluded from the analysis. Input and output matrices were log-transformed, then mean-centered and normalized by standard deviations (column-wise), before performing linear regression [11]. The regression coefficients obtained by this method indicate which parameters the model is most sensitive to with regards to alternans threshold CL, under assumptions of linearity. Linear regression was performed using MATLAB's LinearModel.fit function. Each parameter coefficient was considered significant if the p-value of its t-statistic was greater than 0.05. To evaluate the predictive ability of the regression analysis, we multiplied the regression coefficients by parameter scaling values for the cAF model (log-transformed, mean-centered, and normalized) to obtain the predicted contribution of each parameter to changes in alternans CL [12].