A model for cooperative gating of L-type Ca2+ channels and its effects on cardiac alternans dynamics

In ventricular myocytes, membrane depolarization during the action potential (AP) causes synchronous activation of multiple L-type CaV1.2 channels (LTCCs), which trigger the release of calcium (Ca2+) from the sarcoplasmic reticulum (SR). This results in an increase in intracellular Ca2+ (Cai) that initiates contraction. During pulsus alternans, cardiac contraction is unstable, going from weak to strong in successive beats despite a constant heart rate. These cardiac alternans can be caused by the instability of membrane potential (Vm) due to steep AP duration (APD) restitution (Vm-driven alternans), instability of Cai cycling (Ca2+-driven alternans), or both, and may be modulated by functional coupling between clustered CaV1.2 (e.g. cooperative gating). Here, mathematical analysis and computational models were used to determine how changes in the strength of cooperative gating between LTCCs may impact membrane voltage and intracellular Ca2+ dynamics in the heart. We found that increasing the degree of coupling between LTCCs increases the amplitude of Ca2+ currents (ICaL) and prolongs AP duration (APD). Increased AP duration is known to promote cardiac alternans, a potentially arrhythmogenic substrate. In addition, our analysis shows that increasing the strength of cooperative activation of LTCCs makes the coupling of Ca2+ on the membrane voltage (Cai→Vm coupling) more positive and destabilizes the Vm-Cai dynamics for Vm-driven alternans and Cai-driven alternans, but not for quasiperiodic oscillation. These results suggest that cooperative gating of LTCCs may have a major impact on cardiac excitation-contraction coupling, not only by prolonging APD, but also by altering Cai→Vm coupling and potentially promoting cardiac arrhythmias.

Introduction L-type Ca V 1.2 channels (LTCC) play a critical role in triggering cardiac muscle contraction during the action potential (i.e., excitation-contraction (EC) coupling) [1]. LTCCs are distributed in small clusters of about 10-12 [2][3][4][5] channels along the sarcolemma of these cells [1]. At the membrane potential reached during the plateau phase of the ventricular action potential (AP), LTCCs open, allowing Ca 2+ ions to enter the cell. This Ca 2+ signal is amplified via Ca 2+ -induced Ca 2+ release through opening of ryanodine receptors (RyRs) from the sarcoplasmic reticulum (SR), which causes a cell-wide increase in Ca 2+ that triggers cell contraction [6,7].
Recent experimental studies [2,8,9,28] have suggested that clusters of LTCCs can open and close in unison (i.e., cooperative or coupled gating). Functional coupling between LTCCs requires Ca 2+ for the induction of physical interactions between adjacent channels that ultimately leads to amplification of Ca 2+ influx. This suggests the intriguing hypothesis that cooperative gating of LTCCs may impact membrane voltage (V m ) and intracellular Ca 2+ (Ca i ) cycling dynamics.
In cardiac myocytes, the dynamics of V m are highly nonlinear. The LTCC current (I CaL ) is one of the major currents, which determines the plateau membrane potential and regulates V m dynamics. For example, a slower recovery of LTCC steepens the action potential duration (APD) restitution curve and promotes APD alternans [10][11][12]. Also, reactivation of LTCC during the plateau phase can cause early afterdepolarizations [13][14][15]. Furthermore, intracellular Ca 2+ cycling also has its own nonlinear dynamics [12,16]. In fact, AP clamp experiments showed that the unstable Ca i cycling can induce Ca i transient alternans without APD alternans [16]. The dynamics of Ca i cycling is primary regulated by a steep SR Ca 2+ release vs. SR Ca 2+ load relationship [17][18][19] and thus RyR sensitization or Ca 2+ overload often leads to Ca 2+ alternans.
In this study, we generated a novel computational model of I CaL based on experimental data and APs that incorporates cooperative gating of LTCCs and used it in physiologically detailed AP models to investigate the effects of varied degrees of LTCC coupling on Ca 2+ entry, APD, and the likelihood of promoting voltage and Ca 2+ alternans.

Materials and methods
We developed two computational models. One is a totally stochastic model of cooperative gating (hereafter, we call it 'stochastic model'). In this model, we simulated individual LTCCs described by a stochastic Markov model. Using this model, we investigated I CaL properties by comparing with experimental observations. The other is a deterministic model based on the stochastic model (hereafter, we call it 'ionic model'). Since simulation of the stochastic model is highly computationally intensive and thus time-consuming, we developed the ionic model to permit the investigation of steady state alternans within a practical time frame.

Stochastic model
We used a physiologically detailed subcellular Ca 2+ cycling model as in our previous studies [20][21][22], which is based on the rabbit ventricular myocyte model by Restrepo et al. [23]. Whereas the main structure of our model is similar to the Restrepo model, a key difference in our model is the open probability of RyRs. We reduced the number of RyRs opening during a spark from nearly 100 in the Restrepo model to only 5-10 to fit the experimental observation [20].

Cell geometry
The dimension of this model is 121 μm × 25 μm × 11 μm and there are 19,305 (65×27×11) Ca 2+ release units (CRUs). The separations of CRU are 1.84 μm and 0.9 μm in the longitudinal direction and transverse direction, respectively. Experimental observations indicate that each CRU contains about 10 LTCCs [24,25]. Each CRU in this model contains at least 1 LTCC, but no more than 25. The cluster size obeys Gaussian distribution with a mean of 10 and standard deviation of 3. One where v i is the local cytosolic volume, v s is the local submembrane space volume, v p is the local proximal space volume, v JSR is the local JSR volume, β i is the instantaneous buffer function for c i , β s is the instantaneous buffer function for c s , β p is the instantaneous buffer function for c p , β JSR is the instantaneous buffer function for c JSR , I TCi is time-dependent buffering to Troponin C for c i , I TCs is time-dependent buffering to Troponin C for c s , I dsi is the diffusive current between c s and c i , I up is the uptake current, I ci is the nearest-neighbor diffusive current for c i , I cs is the nearest-neighbor diffusive current for c s , I dps is the diffusive current between c p and c s , I NCX is the sodium-calcium exchange current, I CaBk is the background sarcolemmal membrane Ca flux, I SLCaP is the sarcolemmal membrane Ca pump, I r is the release current, I CaL is the Ltype Ca current, I tr is the JSR refilling current, I cNSR is the nearest-neighbor diffusive current for c NSR , superscript n shows the n-th compartment.
The diffusive currents between different compartments are the same as those previously employed by Restrepo et al. [23].

Cooperative gating of LTCCs and L-type Ca current
As stated above, each CRU contains between 1 and 25 LTCCs. I CaL is described by where i CaL is the single channel current, N L is the number of open channels from 0 to 25, z = VF/(RT). LTCC activity is described by a Markov model with stochastic openings (Fig 1).  In silico analysis of CaV1.2 channel coupling in heart po x and cp x are parameters, which control the sensitivity of the coupling, w 1 and w 2 are the coupling strength. Other rates in Fig 1 are given by: Note that the notations employed here are the same as those used by Mahajan et al. in [26]. A table of constants and the detailed formulation for other current fluxes to reproduce our results can be found in the supplemental material section.

Ionic model
In order to investigate physiological and dynamic regulation of alternans by cooperative gating of LTCCs, we use the action potential (AP) and Ca i cycling model of the ventricular myocyte by Shiferaw et al. [12] since we know which parameters control V m and Ca 2+ dynamics in this model.
The dynamics of membrane voltage (V m ) are described by the equation: where I ion is the total membrane current density, I stim is the stimulus current, and where C m is the cell membrane capacitance. The total membrane current is given by: where I Na is the fast sodium current, I to is the transient outward potassium current, I Kr is the rapid component of the delayed rectifier potassium current, I Ks is the slow component of the delayed rectifier potassium current, I Kp is the plateau potassium current, I K1 is the inward rectifier potassium current, I NaCa is the sodium-Ca 2+ exchanger, and I CaL is the L-type Ca 2+ current. Ca 2+ cycling was modeled by following Shiferaw et al. [19]. This model describes Ca 2+ released from the SR as a summation of local release fluxes distributed throughout the cell. The equations for Ca 2+ cycling are: where c s , c i and c j are the average concentrations of free Ca 2+ in a thin layer just below the cell membrane, in the cytosol, and the SR, with volumes v s , v i and v sr respectively. Here the SR volume includes both JSR and NSR. Also c 0 j is the average JSR Ca 2+ concentration within dyadic junctions in the whole cell. The factors β i and β s describe instantaneous buffering to Calmodulin, the SR membrane, and Troponin C.
All Ca fluxes are divided by v i and have units of μM/ms, which can be converted to units of μA/μF using the conversion factor nFv i /C m , where n is the ionic charge of the charge carrier, and where F is Faraday's constant. Therefore, ionic fluxes can be converted to currents by: where α = Fv i /C m , and where the ion currents are in units of μA/μF.

The L-type Ca 2+ current flux (J CaL ):
J CaL is given by where g Ca is the maximum conductance of J CaL , d is activation, f is voltage-dependent inactivation and f Ca is the Ca 2+ dependent inactivation, i Ca is the single channel current.
To replicate the Markov model shown above, cooperative gating and its Ca 2+ dependence were incorporated in the LTCC model as follows: where γ d is the coupling, which depends on the open probability of LTCC , po x and cs x are parameters, which control the sensitivity of the coupling, w is the coupling strength.
Voltage-dependent inactivation and the Ca 2+ dependent inactivation are given by The single channel current is given by We varied the recovery time constant (τ f ) of the inactivation gate (f) of LTCCs to control the stability of the V m system since it is known that the steepness of APD restitution is sensitive to it. In order to control the stability of the Ca 2+ system, we varied the steepness of the slope (u) of the SR Ca 2+ release function, which controls the sensitivity of release to SR load. The degree of Ca 2+ dependent inactivation (γ) was varied to obtain positive Ca i !V m coupling (γ = 0.7) and negative Ca i !V m coupling (γ = 1.5).
Tables of constants and the detailed formulation for other current fluxes and buffers can be found in the supplemental material section. All programs were written in C++ and run on a 24-node High-Performance Computing cluster and Amazon Cloud Computing Services. An expanded section with all equations can be found in the supplemental material section. The C ++ code is available via our website.

Model of cooperative gating (Stochastic model)
We built a stochastic model of cooperative gating of LTCCs and incorporated this gating modality into the subcellular Ca 2+  Validation of the model involved the use of experimental data. When LTCCs are coupled, simultaneous opening events occurred more often (Fig 2A). Without cooperative gating, events of simultaneous opening of >2 channels are rare. On the other hand, with cooperative gating, events of simultaneous opening of 2 to 6 channels often occurred. Also, open dwell time of the LTCC cluster becomes longer with cooperative gating of LTCCs (Fig 2B vs 2C). We also measured the current-voltage relationship of I CaL and the activation curve. To measure these curves, we used the same protocol used in the experimental study by Dixon et al. [9]. To be more specific, the membrane potential was depolarized from a holding potential of -80 mV to a specified test potential. Fig 3A shows one example of I CaL vs time when V m is depolarized from -80 mV to +20 mV. When cooperative gating is introduced, the peak of I CaL was about 1.5 times larger than that of I CaL without coupling. Cooperative gating of LTCCs shifted the activation curve to the left about 5 or 6 mV ( Fig 3B) and nearly doubled the peak I CaL (Fig 3C). Activation occurs at slightly lower voltage in the model. This discrepancy could be due to species differences as the model is built based on rabbit experiments [26] yet the experimental data collected by Dixon et al. was obtained from mouse cardiomyocytes [9]. Regardless of this, our in silico results are generally consistent with the previously reported experimental observations [9,27]. Furthermore, they support the use of this model to examine the effects of LTCC coupling on Ca 2+ entry, APD, and voltage and Ca 2+ alternans.
To test if cooperative gating promotes alternans, the cell was paced with and without cooperative gating at fast rates. Fig 4 shows the development of alternans and its steady states. The cell was paced at pacing cycle length (PCL) = 300 ms. Without cooperative gating, alternans was not observed (black traces in each panel). However, when cooperative gating was introduced, alternans was developed within 100 beats (Fig 4A: voltage and Fig 4B: Ca 2+ ). It reached the steady state after~30 beats (Fig 4C & 4D). APD alternans amplitude (ΔAPD) is defined as DAPD ¼ ðÀ 1Þ n ðAPD nþ1 À APD n Þ:

Ca 2+ transient alternans amplitude (ΔCa 2+ ) is defined as
DCa 2þ ¼ ðÀ 1Þ n ð½Ca peak nþ1 þ Ca peak n Þ: Without cooperative gating, alternans amplitudes fluctuate around zero. On the other hand, if cooperative gating is introduced, alternans amplitudes stay and fluctuate around certain values (Fig 4E & 4F). When the cell is paced at a faster rate (PCL = 290 ms), alternans occurred in both cases. However, APD and Ca 2+ transient alternans amplitudes were much larger when LTCCs were coupled (S1 Fig). These results are consistent with experimental results [9]. Ca 2+ transient alternans was observed when the cell was paced at PCL = 300 ms using a clamped AP waveform. This implies that Ca 2+ cycling is unstable and contributes development of alternans.

Functional effects of cooperative gating on alternans: Mathematical analysis
Alternans can be caused by instability of V m due to steep APD restitution or instability of Ca i cycling due to steep SR Ca 2+ release vs. SR Ca 2+ load relationship, or both [11,12] (Fig 5A). The V m dynamics and Ca 2+ dynamics are coupled via Ca 2+ -sensitive currents such as I CaL and the Na + -Ca 2+ exchanger (NCX). When the Ca i transient becomes larger, NCX prolongs the APD, whereas the I CaL shortens the APD due to Ca 2+ -induced inactivation of the channel. Therefore, if NCX dominates, APD becomes longer as the Ca i transient becomes larger. In our In silico analysis of CaV1.2 channel coupling in heart previous study [12], we defined this as positive coupling of Ca 2+ on V m (positive Ca i !V m coupling, Fig 5B left). On the other hand, if I CaL dominates, APD becomes shorter as the Ca 2+ transient becomes larger. We defined this as negative coupling of Ca 2+ on V m (negative Ca i !V m coupling, Fig 5B right). The stability of the coupled system is determined by the eigenvalues of a two dimensional map (see ref. [12] for details). The eigenvalues are: where λ v is the eigenvalue associated with the map of the voltage system and λ c is the eigenvalue associated with the map of the Ca 2+ system, and C is coupling of Ca 2+ on the membrane voltage. If V m and Ca i are uncoupled (C = 0), the system becomes unstable (i.e. alternans occurs) when λ v or λ c exceeds unity in absolute value. Black lines in Fig 6A and 6B show stability boundaries (|λ| = 1) when the Ca i !V m coupling is positive (C = 0.1, Fig 6A) and negative (C = -0.1, Fig 6B) from our previous study [12]. In silico analysis of CaV1.2 channel coupling in heart It is also known that Ca 2+ is required for the process of cooperative gating of LTCCs [2,29]. Therefore, larger Ca 2+ transients could tend to prolong APD due to increased strength of cooperative gating of LTCCs. In other words, cooperative gating of LTCCs may promote positive Ca i !V m coupling. Consistent with this, red lines in Fig 6A and 6B show the stability boundary of alternans when Ca i !V m coupling became more positive (C = 0.15, Fig 6A) and less negative (C = -0.05, Fig 6B). An increase in the strength of cooperative gating of LTCCs prolongs APD, which promotes steep APD restitution and increases Ca 2+ influx, which may destabilize Ca 2+ cycling. In addition, this analysis suggests that cooperative gating of LTCCs may destabilize the system (the stable areas became smaller) and promotes alternans except for the quasiperiodic regime.

Functional effects of cooperative gating on alternans: Simulation (ionic model)
To test our theoretical predictions above, we simulated alternans with the ionic model as described in the Methods section. Fig 7A shows that as the coupling strength became larger, APDs became longer. In this simulation, the cell was paced until it reaches the steady state. We chose parameters that do not cause alternans (τ f = 45 ms, u = 3 s -1 ) without coupling (w = 0). Subsequently, the coupling strength was varied and the APD prolongation was measured. The inset shows the action potentials with w = 0.3 and 1.0. When the coupling strength is increased In silico analysis of CaV1.2 channel coupling in heart to 1.0, The peak of the current-voltage curve of I CaL was almost doubled (Fig 7B). Increased coupling strength of LTCCs also shifts the activation curve to the left (Fig 7C). In Fig 7D, I CaL vs time is shown when V m is depolarized from -80 mV to 20 mV. Cooperative gating resulted in a 1.5-fold increase in peak I CaL . These data are consistent with the results with the stochastic model (Fig 3) and experimental observations [9,27].
We also measured the stability boundary to test if increased strength of cooperative gating of LTCCs alters the stability as we predicted in the mathematical analysis (Fig 6). To investigate the effect of the change in channel cooperativity on alternans, we perturbed the cell by changing the coupling strength from w = 0 to w = 0.03. Fig 8A shows the stability boundary when Ca i !V m coupling is positive and Fig 8B shows the stability boundary when Ca i !V m coupling is negative. In both cases, V m -driven alternans and Ca 2+ -driven alternans are promoted with increased coupling strength of LTCCs. On the other hand, quasiperiodic oscillations are not affected. These results are consistent with the theoretical prediction (Fig 6A vs Fig 8A & Fig 6B  vs Fig 8B). We note that this is not due to APD prolongation but due to Ca i !V m coupling. In fact, unlike Fig 8B, simple prolongation of APD by the reduced potassium current stabilizes Ca 2+ -driven alternans when the coupling is negative (S3 Fig). If the coupling strength is within

Cooperative gating, excitation-contraction coupling, and arrhythmias
Cooperative gating of LTCCs facilitates synchronized opening of LTCCs, which may have a major impact on cardiac excitation-contraction coupling due to Ca 2+ signal amplification. In this study, we built stochastic and deterministic computational models of cooperative gating of LTCCs and investigated how this gating modality may affect dynamics of the V m and Ca i cycling system, especially focusing on alternans, which is the arrhythmogenic substrate.  In silico analysis of CaV1.2 channel coupling in heart The novelty of our work is three-fold. Firstly, we model cooperative gating of LTCCs for the first time and add the complexity of this gating phenomenon to the existing models, bringing it more in-line with current thinking on Ca 2+ signaling. We have thus generated a computational model encompassing cooperative gating of LTCCs, which has not been done before. Secondly, we model the effects of cooperative gating on alternans, finding that, in agreement with previously published experimental data [9], aberrant levels of cooperative gating can lead to Ca 2+ alternans. Our theoretical and computational approaches suggest that increases in the strength of cooperative gating of LTCCs promotes positive Ca i !V m coupling and thus promotes V m -driven and Ca-driven alternans. Finally, we confirmed that our model could reproduce experimental data, by specifically examining the effects of changes in the strength of cooperative gating of LTCCs on L-type Ca 2+ currents. The degree of LTCC cooperativity can vary depending on physiological and pathological conditions. Our model provides an in silico means to explore the effects of LTCC cooperative gating under various conditions. In addition, Ca i !V m coupling at the cellular level has been linked to mechanisms of spatially discordant alternans in tissue [30][31][32]. These findings underscore the importance of cooperative gating of LTCCs in excitation-contraction coupling and cardiac arrhythmias.