A bilobal model of Ca²⁺-dependent inactivation to probe the physiology of L-type Ca²⁺ channels

L-type voltage-gated calcium channels (LTCCs) are critical conduits for calcium (Ca²⁺) entry into many excitable cells, including skeletal muscle, cardiomyocytes, and neurons. In ventricular cardiomyocytes, the principal LTCC is the Ca_V1.2 channel, which plays a major role in both shaping the action potential and triggering Ca²⁺ release from the SR, leading to myocyte contraction (Bodi et al., 2005). Disruption of the normal function of this channel has been linked to multiple disorders, including severe forms of long-QT syndrome (LQTS; Crotti et al., 2013; Limpitikul et al., 2014, 2017; Yin et al., 2014; Dick et al., 2016) and heart failure (Ikeda et al., 2009). To properly perform these critical roles, Ca_V1.2 channels use two major forms of feedback regulation: voltage-dependent inactivation (VDI) and Ca²⁺-dependent inactivation (CDI). Of these, CDI is believed to play a dominant role in sculpting action potential morphology (Sun et al., 1997; Morotti et al., 2012). Thus, accurate knowledge of the quantitative details of CDI progression is of critical importance to understanding normal cardiac physiology.

The mediator of CDI in Ca_V1.2 is the ubiquitously expressed Ca²⁺ sensor calmodulin (CaM), which is known to preassociate with the carboxy tail of the channel in its Ca²⁺-free form (apoCaM; Erickson et al., 2001; Pitt et al., 2001; Mori et al., 2004). A channel bound to apoCaM exhibits relatively high open probability upon voltage-mediated activation, known as mode 1 gating (Imredy and Yue, 1994). Upon Ca²⁺ binding, Ca²⁺/CaM departs from its preassociation site, resulting in a reduction of the channel open probability (mode Ca gating; Imredy and Yue, 1994; Zühlke et al., 1999). Ca²⁺/CaM remains associated with the channel via Ca²⁺/CaM effector sites within the amino and carboxy termini of the channel (Ben Johny et al., 2013). Moreover, CaM has the remarkable ability to respond differentially to Ca²⁺ originating from distinct spatial sources. This feature is enabled by the bilobal design of CaM, where a single molecule consists of two semi-independent lobes, each capable of binding Ca²⁺ with different affinities and using differing kinetic schemes to achieve selectivity for a particular source of Ca²⁺ (Dick et al., 2008; Tadross et al., 2008). Specifically, the C-lobe of CaM invariably senses the fluctuation of the local Ca²⁺ concentration at the mouth of the channel, whereas the N-lobe is capable of functioning in two different modes depending on the channel context. For many Ca²⁺ channels, the N-lobe of CaM responds to the cytosolic Ca²⁺ level arising from the cumulative effect of multiple Ca²⁺ sources. However, for Ca_V1.2, the N-lobe acts primarily as a local sensor; although the underlying mechanism for this local selectivity differs from that of the C-lobe (Dick et al., 2008; Tadross et al., 2008). This bilobal scheme plays an important role in physiology. However, the extent of cooperativity between the two lobes remains to be determined.

Disruption of the function of either lobe of CaM has been linked to severe cardiac phenotypes including LQTS and catecholaminergic polymorphic ventricular tachycardia (Crotti et al., 2013; Yin et al., 2014), known as calmodulinopathies. In particular, multiple mutations within the C-lobe of CaM are associated with a frequently lethal form of LQTS (Limpitikul et al., 2014, 2017). In each case, the mutation results in a decrease in Ca²⁺ binding affinity to the C-lobe of CaM, causing a significant reduction in CDI within myocytes harboring the mutation. Thus, full understanding of both normality and pathophysiology requires consideration of the bilobal regulation of CaM on Ca_V1.2.

Given the importance of CDI to physiological function, numerous in silico models have been used to gain a more in-depth understanding of how CDI plays a role in the context of the cardiac action potential. To date, several methods have been deployed to model the CDI of cardiac Ca_V1.2 channels. One method involves treating both VDI and CDI as canonical Hodgkin and Huxley gates, dependent on voltage or Ca²⁺ concentration, respectively; that is, total Ca²⁺ current is a product of conductance, VDI “gate,” CDI “gate,” and a driving force based on the Goldmann–Hodgkin–Katz equation (Hund and Rudy, 2004; ten Tusscher et al., 2004; ten Tusscher and Panfilov, 2006). Building on this, a second method includes Markovian models containing closed, open, and inactivated states. Here, rate constants for the transition into a fraction of the inactivated states depend on voltage (representing VDI), whereas transitions into the remaining inactivated states depend on Ca²⁺ concentration (representing CDI). Importantly, most of these models do not incorporate the allosteric mechanism of CDI, where Ca²⁺-inactivated channels maintain the ability to open, but with a reduced (yet nonzero) open probability (Faber et al., 2007; Mahajan et al., 2008). A few models have partially addressed this issue; however, even in most of these models, the majority of inactivation is carried by VDI, in contrast to the known physiological mechanism of Ca_V1.2 channels (Jafri et al., 1998; Decker et al., 2009). In attempting to capture these distinct CDI kinetics, some have included the addition of two time constants for CDI (Hund and Rudy, 2004; O’Hara et al., 2011), yet the bilobal response of Ca_V1.2 remains unexplored. Moreover, Greenstein and Winslow (2002) and Greenstein et al. (2006) demonstrated that the proper physiological balance between CDI and VDI (in which LTCC inactivation is strongly dominated by CDI), and hence the proper relationship between CDI strength and action potential morphology, can be achieved in simulations in which LTCCs are localized to independently functioning cardiac dyads within a myocyte. Although all of these approximations of Ca_V1.2 inactivation can indeed produce Ca²⁺ current (I_Ca) with kinetics grossly resembling experimental I_Ca recorded during an action potential, they do not take into account the allosteric and bilobal nature of CaM regulation of Ca_V1.2, thus limiting their utility in describing Ca_V1.2 gating in response to variable spatial calcium signals (Sun et al., 1997; Dick et al., 2008; Tadross et al., 2008).

One obstacle to the development of accurate CDI models is a paucity of experimental data. As Ca_V1.2 becomes permeable to the Ca²⁺ ions that modulate its own opening, it is challenging to accurately control and measure the amount of Ca²⁺ that preassociated CaM actually senses. Moreover, during an action potential, there is at least a ∼10-fold difference in Ca²⁺ concentration between the dyadic space facing the channel pore and the rest of the cytoplasm (Neher, 1998). Most experimental values of Ca²⁺ used to drive or validate the CDI models are measured by cytosolic Ca²⁺-sensitive fluorescence dyes, where the readout represents the average of bulk (global) and subspace (local) Ca²⁺ concentrations. To accurately measure the Ca²⁺ sensed by a CaM residing on a channel, we take advantage of photo-uncaging (Lipp et al., 1996) to deliver a controlled amount of Ca²⁺ to the channels, while simultaneously monitoring the cytosolic concentration of Ca²⁺, a method that has been validated in the context of both L-type and P/Q-type Ca²⁺ channels (Tadross et al., 2013; Lee et al., 2015; Dick et al., 2016). As we allow the Ca²⁺ cage to fully dialyze before the photolysis, this Ca²⁺ measurement will accurately represent both the bulk and local Ca²⁺ concentration.

Using this technique, we recorded CDI for each lobe of CaM individually and resolved the kinetic response to measured levels of Ca²⁺. This new high-resolution data enabled the creation of a four-state model for independent N- or C-lobe–mediated CDI. Our models successfully recapitulate both the kinetic and steady-state properties of CDI for each lobe. Next, we considered the integrated effect of bilobal CaM on Ca_V1.2 CDI. Using experimental CDI measurements for wild-type CaM (CaM_WT), we constructed a 16-state model of bilobal CaM regulation. Not only can this full model effectively capture both the kinetics and steady-state behavior of CDI under physiological conditions, but it also recapitulates the behavior of calmodulinopathic CaM. Thus, this CDI model could prove to be an invaluable asset in understanding both physiological and pathological states of Ca_V1.2 regulation.

The Ca_V1.2 channel used in this study was made up of the human α_1C cardiac isoform corresponding to Acc# Z34810 within the pcDNA3 vector. The original channel backbone was a gift from T.W. Soong (National University of Singapore, Singapore; Tang et al., 2004) and was modified as previously described (Dick et al., 2016) to conform to the desired cardiac isoform. The point mutation E736A (in the S6 helix of domain II; Ellinor et al., 1995) was introduced using standard overlap extension PCR, resulting in a Ca²⁺ channel with increased selectivity for Li⁺ (Ca_V1.2^E2A; Ellinor et al., 1995; Dick et al., 2016). The altered selectivity and decrease of pore block for this channel was previously determined (Ellinor et al., 1995; Dick et al., 2016).

Engineered rat brain CaM₁₂, CaM₃₄, and CaM₁₂₃₄, in which the aspartate-to-alanine mutations were introduced in the EF hands of N-lobe, C-lobe, and both lobes, respectively, were gifts from John P. Adelman (Vollum Institute, Portland, OR). The LQTS CaM_D130G mutation was generated using QuikChange site-directed mutagenesis (Agilent) on rat brain CaM (M17069) in the pcDNA3 vector (Invitrogen). Both WT and engineered CaMs were cloned into the pIRES2-ECFP vector, modified from pIRES2-EGFP (Clontech Laboratories), using NheI and BglII, thus allowing visualization of positively transfected cells.

HEK293 cells were cultured on glass coverslips in 10-cm dishes. α_1C along with its auxiliary subunits β_2a (M80545) and α₂δ (NM012919.2) were transiently transfected along with either WT or mutant CaM using a standard calcium phosphate method (Peterson et al., 1999). 8 µg of each plasmid was used, along with 2 µg of simian virus 40 T antigen complementary DNA to ensure robust expression.

Whole-cell voltage-clamp recordings of HEK293 cells were performed at room temperature 1–2 d after transfection. Recordings were obtained using an Axopatch 200B amplifier (Axon Instruments). Whole-cell voltage-clamp records were digitally sampled at 10 kHz and low-pass filtered at 5 kHz. P/8 leak subtraction was used, with series resistances of 1–2 MΩ. Experiments were performed in external solutions containing (in mM) 80 TEA-MeSO₃, 10 HEPES, pH 7.4, 80 LiCl, and 2 EGTA at 300 mOsm adjusted with glucose. Internal solutions contained (in mM) 135 CsCl, 40 HEPES, pH 7.4, 1–2 citrate, either 5 µM Fluo-2 (TefLabs) + 5 µM Fluo-2 Low Affinity (TefLabs) or 10 µM Fluo-4FF (Invitrogen), 2.5 µM Alexa Fluor 568 (Invitrogen), 2 DMNP-EDTA (Invitrogen), and 0.25–1 CaCl₂. Measurement of CDI was done in response to a 0-mV step depolarization, which is close to the peak of the current–voltage relationship for all combinations of channel and CaM coexpression (Supplementary text, Co-expression of engineered CaM).

During current recordings, Ca²⁺ concentration was simultaneously measured optically using Ca²⁺ indicator dyes. To ensure proper quantification of Ca²⁺ at various levels, we selected Ca²⁺ dyes with variable Ca²⁺ affinities as required for the range of Ca²⁺ measured. For all Ca²⁺ steps <10 µM, we used an internal solution containing the mixture (see above) of Fluo-2 Low Affinity (K_d = 6.7 µM) and Fluo-2 High Affinity (K_d = 230 nM). This allowed clear resolution of low basal Ca²⁺ (pre-flash) as well as the ability to distinguish even small or moderate Ca²⁺ steps (post-flash; Lee et al., 2015). This enhanced range of Ca²⁺ measurements permitted precise measurement of the time-varying Ca²⁺ concentration after the flash, allowing our model to account for the kinetic changes in both Ca²⁺ and CDI. To obtain accurate measurements of Ca²⁺ steps in the range >10 µM, Fluo-4FF (K_d = 9.7 µM) was used (Dick et al., 2016). To decrease the variability across experiments, each dye mixture was premixed in a large stock and used throughout the study. Calibration of premixed stocks of Fluo-2/Alexa was performed by Lee et al. (2015) and Fluo-4FF/Alexa by Dick et al. (2016). The dyes were excited with a 514-nm argon laser, via a 545DCLP dichroic mirror and either a 545/40BP (Fluo) or 580LP (Alexa Fluor 568) filter. Fluorescence intensity was measured via photomultiplier tubes at a sampling frequency of 10 kHz. Brief UV pulses used to produce Ca²⁺ steps were generated by a Cairn UV flash photolysis system. The amount of Ca²⁺ delivered to channels was controlled by adding different ratios of CaCl₂ to DM-nitrophen (DMN) in the pipette solution and varying the UV intensities achieved by adjusting the voltage (50–200 V) and total capacitance (500–4,000 μF) of the flash photolysis system.

Data are presented as mean ± SEM. Statistical analysis of the data was performed using Prism 5 (GraphPad). To compare the CDI response to the level of intracellular Ca²⁺ concentration for each experimental condition, CDI as a function of Ca²⁺ was fitted with Hill equations, and the extra sum-of-squares F test was used to determine statistically significant difference among the K_d values.

Experimentally measured or apparent CDI was a combination of true CDI and residual Ca²⁺ pore block of Ca_V1.2^E2A channel (Fig. 1, A–D). Careful characterization of this residual pore block allows quantification of the magnitude of true CDI (Supplementary text, Derivation of apparent and true CDI) which is used in all computational simulations. For N-lobe, C-lobe, and bilobal models of CDI, a set of ordinary differential equations (ODEs) was used to represent the probability of being in each state (Figs. 3 A, 4 A, and 5 A), with the sum of all probabilities constrained to be 1. Detailed equations are described in the Supplementary text, Differential equations for N-lobe and C-lobe CDI and Differential equations for bilobal CDI. All simulations were performed using Matlab (MathWorks). ODEs were solved using Matlab function ode23t. The measured Ca²⁺ concentration was fitted with a double exponential equation, and this idealized Ca²⁺ signal was used as a model input (Supplementary text, Model setup and parameter optimization). Model optimization was performed by using Matlab function patternsearch to minimize the error between experimental and simulated current (Supplementary text, Model setup and parameter optimization). Parameter constraints for all models are described in detail in Supplementary information, Parameter constraints.

The supplemental text describes the effect of co-expression of engineered CaM on voltage dependence of Ca_V1.2 channel activation; contains derivation of true and apparent CDI; contains differential equations describing N-lobe and C-lobe CDI; contains describing differential equations for bilobal CDI; contains model setup and a detailed method for parameter optimization; describes parameter constraints; demonstrates the importance of cooperativity factors in the bilobal CDI model; contains derivation of cooperativity between CaM lobes; and describes parameter sensitivity. Fig. S1 shows the current–voltage relationship for Ca_V1.2^E2A. Fig. S2 shows a schematic of CDI measurement as a function of [Ca²⁺]. Fig. S3 shows a method for idealization of [Ca²⁺]. Fig. S4 shows a cost function schematic. Fig. S5 shows model results in the absence of cooperativity. Fig. S6 shows evaluation of parameter sensitivity in the N-lobe model. Fig. S7 shows evaluation of parameter sensitivity in the C-lobe model. Fig. S8 shows evaluation of parameter sensitivity in the bilobal model. Table S1 shows [Ca²⁺] parameters used in model simulations. Table S2 shows quantification of the sensitivity of model parameters.

To create a detailed kinetic model of CDI, we first need to obtain high-resolution CDI data for each lobe of CaM at known Ca²⁺ concentrations. To do this, we must overcome the experimental challenge of coupled Ca_V gating and Ca²⁺ entry, which we accomplished by adapting a previously described photo-uncaging method (Tadross et al., 2013; Lee et al., 2015; Dick et al., 2016). We use Li⁺ as the charge carrier across the channel (Fig. 1 A), while controlling intracellular Ca²⁺ via the Ca²⁺ cage DMN. However, because of the high affinity of Ca²⁺ for the channel pore, intracellular Ca²⁺ has been shown to block the entry of monovalent ions in a process known as pore block (Ellinor et al., 1995). To mitigate this artifact, we used a previously described Ca_V1.2 channel in which the mutation E736A was inserted into the selectivity filter of domain II in the α₁ subunit of the Ca_V1.2 channel (Ca_V1.2^E2A), which shifts the selectivity of the channel in favor of Li⁺ (Ellinor et al., 1995; Dick et al., 2016).

We are now able to deliver a spatially uniform concentration of intracellular Ca²⁺ such that measured bulk cytosolic Ca²⁺ is identical to that surrounding the CaM residing on the channel (Tadross et al., 2013; Lee et al., 2015; Dick et al., 2016). Specifically, during patch-clamp recording, we dialyzed DMN into HEK293 cells cotransfected with CaM and human cardiac Ca_V1.2 channels. Then, upon delivery of a UV pulse, the ester bonds of DMN are cleaved and Ca²⁺ is released into the cytosol. In this setting, the magnitude of the Ca²⁺ step delivered can be controlled by adjusting the intensity of the UV light and varying the ratio between Ca²⁺ and DMN in the dialysate (Figs. 1, B and C; Tadross et al., 2013; Lee et al., 2015; Dick et al., 2016). To measure the real-time concentration of Ca²⁺, we added Ca²⁺-sensitive dyes to the dialysate and simultaneously recorded the fluorescence intensity and Li⁺ current (Fig. 1 A). Because the physiologically relevant range of Ca²⁺ is larger than the dynamic range of a single Ca²⁺ dye, we combined known ratios of low- and high-affinity Ca²⁺ dyes (Fluo-2 HighAff and Fluo-2 LowAff) to accurately measure [Ca²⁺] across the entire range (Lee et al., 2015).

To characterize the extent of any residual Ca²⁺ pore block in Ca_V1.2^E2A, we coexpressed the channel with CaM₁₂₃₄, an engineered CaM in which the Ca²⁺ binding of all four EF hands has been disrupted (Dick et al., 2016). Under these conditions, Li⁺ current displays a slow decay in current amplitude, indicative of VDI (Fig. 1 A). Upon Ca²⁺ photo-uncaging, a small instantaneous drop in current demonstrates the remaining pore block of the Ca_V1.2^E2A channel (Fig. 1, B and C). This residual pore block can be quantified as:

where I_preUV is current measured 20 ms before the UV pluses, I_end is current remaining at the end of the voltage step, and nf (no flash) and f (flash) represent current without and with UV flash delivered 150 ms after channel activation, respectively. For Ca_V1.2^E2A, f_block increases monotonically as a function of [Ca²⁺] (Fig. 1 D) and can be fitted with the Hill equation,

with amplitude (amp_block) = 1, Hill coefficient (nHill) = 1, and K_d/block = 33 µM for Fluo-2/Alexa or 55 µM for Fluo-4FF/Alexa, which are comparable to previously described values (Ellinor et al., 1995; Dick et al., 2016).

This setup can now be used to probe the extent of CDI of Ca_V1.2 at various cytosolic levels of Ca²⁺. When Ca_V1.2^E2A is now coexpressed with wild-type CaM (CaM_WT), a step release of Ca²⁺ (Fig. 1 F, red) results in robust inactivation of the Li⁺ current (Fig. 1 F, black) indicating CDI. Importantly, the confounding effects of VDI can be eliminated in two ways. First, we coexpressed our channels with the β_2a subunit, which is known to minimize the extent of VDI in the channel complex (Dafi et al., 2004). Second, we can eliminate any residual VDI by using the preflash current as our baseline for CDI measurements (Eq. 3). Further, by varying the amplitude of the Ca²⁺ step, we can quantify the extent of steady-state CDI as a function of Ca²⁺ (Fig. 1 G). Here, the normalized ratio of current reduction after Ca²⁺ delivery (f_CDI,app, Fig. 1 G) represents the combination of true CDI (f_CDI,true, Fig. 1 H) and residual pore block (f_block, Fig. 1 D):

(see Supplementary text, Derivation of apparent and true CDI). Thus, the true CDI can be calculated (Fig. 1 H) and fitted with a Hill equation of the form

where CDI_max,WT = 0.88, nHill_WT = 2.5, and K_d,WT = 0.52 µM (R² = 0.97), comparable to previously described values (Dick et al., 2016). Thus, we have a robust method for quantifying the extent of CDI in response to controlled concentrations of intracellular Ca²⁺, providing data with which to optimize our in silico models.

Prior studies have demonstrated separate forms of CDI mediated independently by the N- or C-lobe of CaM. Importantly, these two forms of CDI can respond to spatially distinct Ca²⁺ sources, each using different mechanistic schemes (Dick et al., 2008; Tadross et al., 2008). However, this unique feature of CaM is not represented in published in silico models of CDI. This lack of lobe-specific detail is, in large part, the result of limited data pertaining to each lobe of CaM. We therefore undertook detailed experimental characterization of N- and C-lobe–mediated CDI.

To begin, we use an engineered CaM (Keen et al., 1999), in which Ca²⁺ binding of the C-lobe is abolished (CaM₃₄), allowing the resolution of N-lobe–mediated CDI alone (Peterson et al., 1999). Ca_V1.2^E2A coexpressed with CaM₃₄ displays robust Li⁺ current at baseline (Fig. 2 A). However, upon photo-uncaging of Ca²⁺, the current decays exponentially (Fig. 2 B, black), with notably different kinetics compared with CaM_WT (Fig. 1 F). By varying the amount of Ca²⁺ delivered, we can describe the extent of steady-state N-lobe–mediated CDI as a function of intracellular Ca²⁺ (Fig. 2 C). This relationship can be fitted with a Hill equation (Eq. 4) with CDI_max,N = 0.62, Hill coefficient nHill_N = 1.8, and K_d,N = 0.9 µM (R² = 0.79).

Similarly, CDI mediated by the C-lobe CaM can be probed by measuring Li⁺ current through Ca_V1.2^E2A coexpressed with an engineered CaM with restricted C-lobe–only Ca²⁺ binding (CaM₁₂; Keen et al., 1999). Here, photo-uncaging Ca²⁺ results in a rapid inactivation of the Li⁺ current (Fig. 2 E, black) at a faster rate than in the presence of CaM₃₄ (Fig. 2 B). Notably, this increased kinetic response of C-lobe CaM occurs at each level of Ca²⁺ measured. Moreover, steady-state CDI mediated by the C-lobe of CaM can also be fitted by a Hill equation (Eq. 4), with CDI_max,C = 0.66, Hill coefficient nHill_C = 1.4, and K_d,C = 1.15 µM (R² = 0.93; P > 0.05 compared with K_d,N). Importantly, the K_d values obtained for N- and C-lobe CaM are each significantly larger than that obtained for CaM_WT (0.9, 1.15, and 0.5 µM respectively, P < 0.0001).

Having achieved robust experimental measurements of CDI mediated by each lobe of CaM, we now have the necessary components to build a kinetic model of CDI in silico. We begin by first considering each lobe of CaM independently, using the experimental results obtained for CaM₃₄ and CaM₁₂ to optimize the model parameters. Because current LTCC models sufficiently address voltage-dependent activation and the voltage-dependent component of inactivation (VDI), we restrict our model to the context of CDI. Importantly, this focus does not ignore the importance of VDI, as this component can be added into any model owing to the independence of these two processes (Luo and Rudy, 1994; Jafri et al., 1998; Hirano and Hiraoka, 2003; Findlay et al., 2008), as such:

where I(t) is the total current remaining after CDI and VDI are considered, and I_max(t) is the maximum current amplitude based solely on channel activation and permeation. Thus, our CDI models are formulated to describe the time course of CDI after channels have fully activated.

To create kinetic models of CDI, we first considered the allosteric nature of this process (Faber et al., 2007; Mahajan et al., 2008). That is, fully activated channels bound to apoCaM open with a relatively high P_o (mode 1 gating), whereas channels harboring a Ca²⁺-bound CaM transition into a low P_o mode of gating “mode Ca” (Imredy and Yue, 1994; Adams et al., 2014). Once in mode Ca, channels maintain the ability to open, but with a reduced open probability, resulting in CDI. Importantly, recent work indicates that the transition to mode Ca occurs upon the departure of apoCaM, such that channels devoid of apoCaM may be indistinguishable from channels bound to Ca²⁺/CaM (Adams et al., 2014). We first applied these principles to create a model of N-lobe–mediated CDI. Here, we assume that the C-lobe will remain in the apo-state and thus will not contribute to CDI, whereas the N-lobe will independently orchestrate its own form of CDI (Fig. 2 C). We therefore created a four-state model of N-lobe CDI, such that each individual state represents a channel in either mode 1 (apoCaM bound, high P_o) or mode Ca (empty channel or Ca²⁺/CaM bound, low P_o; Fig. 3 A). Thus, our model represents the minimal number of states required to reproduce CDI based on currently available structure–function data (Imredy and Yue, 1994; Tadross et al., 2008, 2010; Ben Johny et al., 2013; Adams et al., 2014). Specifically, state 1 represents apoCaM₃₄ bound to the carboxy terminus of the channel. State 2 depicts a departure of the N-lobe from the channel, which then allows Ca²⁺ to bind transitioning to state 3. The calcified N-lobe can then subsequently bind to its effector site on the amino terminus of the channel, represented by state 4 (Dick et al., 2008; Tadross et al., 2008). For the binding of Ca²⁺ to CaM, we assumed simultaneous Ca²⁺ (un)binding, as has been demonstrated for the two EF hands within each lobe of CaM (Martin et al., 1985; Linse et al., 1991). Moreover, we define the open probability of state 1 (apoCaM₃₄ bound) to be P_O,max while assuming the open probabilities of states 2, 3, and 4 are reduced equally by a factor f_N (fractional open probability in mode Ca (Adams et al., 2014), such that and the open probability of each of these mode Ca states is f_N * P_O,max. From our experimental measurements (Fig. 2 C), we obtain f_N = 0.38. We can then express the time-varying evolution of CDI as

where M_x(t) is the probability of channels being in state x at time t (see Supplementary text, Differential equations for N-lobe and C-lobe CDI, for detailed model equations). N-lobe states were defined as shown in Fig. 3 A, subject to the constraints as described in Supplementary text, Parameter constraints. Overall, the model contains six free parameters; however, constraints were placed on k_on,N and k_off,N based on prior experimental reports (Bayley et al., 1984; Martin et al., 1985; Stemmer and Klee, 1994; Black et al., 2005; Saucerman and Bers, 2008), leaving only four free parameters requiring full optimization (Supplementary text, Model setup and parameter optimization, and Parameter constraints).

We next optimized our model parameters (Table 1) to minimize the error between our model output and our kinetic (Fig. 2 B) and steady-state (Fig. 2 C) experimental data for CaM₃₄. To begin, we fitted our experimental Ca²⁺ measurements with a double exponential equation to provide a smooth Ca²⁺ input waveform for our model (Fig. 3 B, bottom; Supplementary text, Model setup and parameter optimization). Accordingly, we simulated currents in response to three different Ca²⁺ step inputs, yielding optimized rate constants as described in Table 1. Although these parameters do not definitively quantify all underlying aspects of the biological system, our model effectively recapitulates the kinetic and steady-state aspects of our experimental data (Figs. 3, B and C; red, simulation; black, experiment), with k_on,N and k_off,N comparable to experimentally measured values (Bayley et al., 1984; Martin et al., 1985; Stemmer and Klee, 1994; Black et al., 2005; Saucerman and Bers, 2008). Moreover, the parameters of this N-lobe–mediated CDI model are consistent with the known mechanisms of spatial Ca²⁺ selectivity of N-lobe CaM (Supplementary text, Parameter constraints; Dick et al., 2008; Tadross et al., 2008).

Having established a functional model of N-lobe CaM signaling, we next sought to create a similar model for C-lobe–mediated CDI. Using the same approach as with the N-lobe model, we created a four-state model of C-lobe mediated CDI (Fig. 4 A, see Supplementary text, Differential equations for N-lobe and C-lobe CDI, for detailed differential equations for each mode). We again defined the open probability of state 1 (apoCaM₁₂ bound) to be P_O,max while assuming the open probability of states 2, 3, and 4 to be f_C * P_O,max, with f_C defined as From our experimental measurements (Fig. 2 F), we obtain f_C = 0.34. We then express the time-varying extent of CDI as

where M_x(t) is the probability of channels being in state x at time t (see Supplementary text, Differential equations for N-lobe and C-lobe CDI, for detailed model equations). We then optimized the free parameters (Table 1) fitting the model to the kinetic (Fig. 2 E) and steady-state (Fig. 2 F) experimental data for CaM₁₂. The simulated currents (Fig. 4 B) reproduce experimental results well, capturing the initial rapidly inactivating phase, the following slowly inactivating phase, and the steady-state response to Ca²⁺ (Fig. 4 C). Again, the kinetics of Ca²⁺ binding and unbinding to C-lobe of CaM, k_on,C and k_off,C, are comparable to those experimentally described for this lobe (Bayley et al., 1984; Martin et al., 1985; Stemmer and Klee, 1994; Black et al., 2005; Saucerman and Bers, 2008). Moreover, the parameters of this C-lobe–mediated CDI model is consistent with the invariably local Ca²⁺ selectivity of C-lobe CaM (See Supplementary text, Parameter constraints, for details; Dick et al., 2008; Tadross et al., 2008).

Having established independent models of lobe-specific CaM-mediated CDI, we next sought to integrate these models into a full model representing the kinetic response of both lobes of CaM simultaneously (CaM_WT). Cognizant of the potential interdependence of the two lobes of CaM, we created a bilobal model of CDI consisting of 16 states which represent all possible combinations of the minimal four-state N- and C-lobe models (Fig. 5 A). Here, horizontal transitions represent C-lobe state changes, and vertical transitions represent N-lobe state changes. Our model incorporates the optimized rate constants from both the N- and C-lobe models (Figs. 3 A and 4 A and Table 1), as well as multiple cooperativity factors (g, h, k, l, m, n, p, and q) to account for potential communication between the lobes (Fig. 5 A). Consideration of thermodynamic microscopic reversibility reduces this to four free parameters (Table 1) as described in Supplementary Information, Differential equations for bilobal CDI and Parameter constraints. In addition, we assumed that a transiently dissociated CaM (states 6, 7, 10, and 11) is unlikely to diffuse away from the channel alcove owing to the limiting rate of exchange between cytosolic and channel-associated CaMs (Chaudhuri et al., 2005). As in the reduced models, the open probability in state 1 (both lobes of apoCaM bound) is defined to be P_O,max; the open probability in states 2, 3, and 4 is defined as f_N * P_O,max; the open probability in states 5, 9, and 13 is f_C * P_O,max; and the open probability in states 6, 7, 8, 10, 11, 12, 14, 15, and 16 is f_WT * P_O,max, where f_WT = 1 − CDI_max = 0.12 (Fig. 1 H) and CDI is defined as: where M_x(t) is the probability of channels being in state x at time t (see detailed model equations in Supplementary text, Differential equations for bilobal CDI).

Fitting this model to our experimental data for Ca_V1.2^E2A coexpressed with CaM_WT enabled us to optimize the parameters (Table 1, Supplementary text, Differential equations for bilobal CDI and Parameter constraints), resulting in a bilobal model that approximates both the kinetic (Fig. 5 B) and steady-state (Fig. 5 C) behavior of CaM-mediated CDI. In addition, the model predicts cooperativity factors with both positive (>1, increased likelihood of the transition) and negative (<1, decreased likelihood of the transition) cooperativity. Importantly, setting all cooperativity factors to unity (no cooperativity) results in a model that fails to reproduce the experimental data (Supplementary text, The importance of cooperativity factors). Therefore, although these cooperativity factors may not quantify the precise physiological mechanism, some form of cooperativity is likely to exist. Thus, our model explicitly captures the mechanistic interplay between the two lobes of CaM, an important feature of CDI not found in previous models.

To gauge the ability of our bilobal CDI model to make predictions under pathophysiological conditions, we applied it in the context of a known calmodulinopathic mutation D130G. This missense mutation occurs within the fourth EF hand of CaM, resulting a ∼50-fold decrease in the Ca²⁺ binding affinity of the C-lobe (Crotti et al., 2013; Hwang et al., 2014) and results in a significantly increased QT interval, rendering patients susceptible to arrhythmias (Arevalo et al., 2007). At the cellular level, this is demonstrated by a prolongation of the action potential (Fig. 6 A), which is largely caused by disruption of Ca_V1.2 CDI (Limpitikul et al., 2014). We thus applied our 16-state model to CaM_D130G-mediated CDI. To begin, we experimentally probed CDI from Ca_V1.2 channels that coexpressed CaM_D130G using our photo-uncaging technique (Fig. 6, B and C, black). Here, the CDI is significantly reduced even at very high concentrations of Ca²⁺ (Fig. 6 C). Moreover, the steady-state CDI mediated by CaM_D130G (Fig. 6 C) can be fitted with a Hill equation (Eq. 4) with CDI_max,D130G = 0.62, nHill_D130G = 1.8, and K_d,D130G = 0.9 µM (R² = 0.98). Notably, this represents the identical parameters used to describe the steady-state behavior of N-lobe–mediated CDI. Thus, this mutation in a single EF hand appears to disrupt CDI to the same extent as CaM₃₄, in which both C-lobe EF hands are ablated.

To simulate this mutation with our 16-state CDI model, we decreased k_on,C by 50-fold to account for the loss of Ca²⁺ binding affinity for this lobe and dramatically reduced the cooperativity factors g and m, preventing CDI from proceeding via C-lobe CaM (Table 2). The resultant model was able to predict the kinetic (Fig. 6 B) and steady-state (Fig. 6 C) properties of CaM_D130G-mediated CDI, thus demonstrating the utility of this bilobal CDI model in both physiological and pathological conditions.

CDI is a critical regulatory feature of Ca_V1.2 channels, playing a pivotal role in setting the action potential duration in the heart. Although multiple cardiac and neuronal models use CDI components (Hund and Rudy, 2004; ten Tusscher et al., 2004; Faber et al., 2007; Mahajan et al., 2008), few take into account the detailed kinetics or bilobal nature of CaM binding to the channel. Here, we obtain rigorous experimental data and use it to create a detailed kinetic model of Ca_V1.2 CDI featuring full bilobal CaM regulation. Specifically, we measured the kinetic response of channel inactivation caused by measurable changes in intracellular Ca²⁺ by pairing photo-uncaging with simultaneous patch clamp recording (Tadross et al., 2013; Lee et al., 2015; Dick et al., 2016). This method allowed us to uncouple channel gating from intracellular Ca²⁺ dynamics, enabling precise control of bulk cytosolic Ca²⁺ and overcoming the ambiguity caused by diffusion-limited Ca²⁺ dispersion in traditional patch-clamp experiments. Moreover, by using single and combinations of Ca²⁺-sensitive dyes with different affinities, we were able to accurately measure [Ca²⁺] across a wide range of values from resting (∼0.1 µM) to locally elevated (as high as 100 µM at the mouth of the channel; Neher, 1998). As such, we concurrently measured CDI in response to known concentrations of Ca²⁺. This use of high-resolution data enabled the development of two independent single-lobe models of Ca_V1.2 CDI and a 16-state bilobal model of CDI.

To build this bilobal model of CDI, we used a reductionist approach, first examining the CDI elicited by each lobe of CaM independently. Our experimental data showed that despite minimal difference in steady-state CDI profiles, N-lobe-mediated CDI exhibits much slower kinetics compared with C-lobe-mediated CDI. Using this data, we created a minimal four-state single-lobe model for each lobe of CaM and incorporated these two models into a 16-state bilobal model of CDI, enabling us to gauge the cooperativity between the two lobes in the context of Ca_V1.2 binding. Although cooperativity between the N- and C-lobe of CaM is well known in response to Ca²⁺ binding (Teleman et al., 1986; Yazawa et al., 1987; Jaren et al., 2002), the extent of cooperativity in the initiation of Ca_V1.2 CDI remains undetermined. Given the importance of this feature, we can corroborate the positive cooperativity factor implied by our model by scrutinizing our experimental data. Specifically, we express the closed-form solution of steady-state bilobal mediated CDI as

where K_N and K_C are the K_d of N- and C-lobe CaM respectively; P_N, P_C, and P_WT are the relative open probabilities of mode Ca compared with mode 1 for N-, C-, and both lobes of CaM, respectively; and λ is the cooperativity factor (see Supplementary text, Derivation of cooperativity between CaM lobes, for derivation). Here, λ >1 and λ < 1 represent positive and negative cooperativity between CaM lobes, respectively. Applying this equation to our steady-state CaM_WT CDI data (Fig. 1 H) reveals λ = ∼5, demonstrating that both lobes of CaM modulate Ca_V1.2 in an interdependent manner with positive cooperativity, which implies that independent four-state N- and C-lobe models might not be sufficient to correctly describe the kinetics of the bilobal CDI process. Moreover, cooperativity implies that perturbation in one lobe of CaM can affect overall CDI to a larger extent than would be predicted in a model in which each lobe acts independently and may be an important feature in understanding the pathogenesis of diseases such as calmodulinopathies.

We therefore applied our bilobal CDI model to simulate the effect of the LQTS-associated calmodulinopathy mutation, D130G (Crotti et al., 2013; Limpitikul et al., 2014, 2017; Yin et al., 2014). Residing in the fourth EF hand of CaM, this mutation minimally increases the apoCaM binding affinity to Ca_V1.2 (approximately two-fold; Limpitikul et al., 2014) while dramatically diminishing the Ca²⁺ binding affinity of C-lobe CaM (∼50-fold; Crotti et al., 2013; Hwang et al., 2014), with no discernable effect on the Ca²⁺ binding affinity of N-lobe CaM (Hwang et al., 2014). Interestingly, although the C-lobe of CaM_D130G maintains some Ca²⁺-binding capacity, our experimental data show that overall, CaM_D130G behaves nearly identically to CaM₃₄ (Figs. 2 C and 6 C). Moreover, prediction from our bilobal CDI model corroborates this finding. Upon decreasing k_on,C by 50-fold and dramatically reducing the cooperativity factors g and m (thus preventing CDI from proceeding via C-lobe CaM), the model accurately recapitulates both the kinetic and steady-state behavior of CaM_D130G. If strong interlobal cooperativity could be verified experimentally, it would imply that when the Ca²⁺-binding affinity of one lobe of CaM falls under a certain threshold, CDI mediated by that lobe may become physiologically ineffective. In fact, experimental data shows that even a fivefold reduction in the C-lobe Ca²⁺ binding affinity (F142L mutation; Crotti et al., 2013; Hwang et al., 2014) can alter the overall behavior of CaM to resemble that of CaM₃₄ (Limpitikul et al., 2014).

Although our models can approximate the kinetics of CDI in both physiological and pathological conditions relatively well, there exist experimental and computational limitations. For example, rate constants for the CDI models were optimized using data obtained at room temperature, which are likely slower than those at a physiological 37°C. However, with known temperature coefficient (Q₁₀), rate constants at 37°C can be extrapolated (Faber et al., 2007; Mahajan et al., 2008). In addition, here, Ca_V1.2 channels and CaM variants were heterologously expressed in HEK293 cells, which lack some modulatory mechanisms inherent to cardiomyocytes, such as adrenergic regulation of Ca_V1.2 (Miriyala et al., 2008). Although we were not able to incorporate all regulatory feedbacks of Ca_V1.2, this reductionist approach could be advantageous in that it allows the CDI module to be incorporated into models of different cell types. In addition, little is known about the conformational changes that result in CDI and thus cannot be explicitly included in the model. Moreover, we used the assumption that states 2 and 3 will have a low open probability based on the lack of apoCaM binding; however, it must be noted that this hypothesis is based on the work of Adams et al. (2014) and remains to be fully verified. Furthermore, consideration of the sensitivity of the model to each parameter (Supplementary text, Parameter sensitivity) reveals that the final parameter set does not comprise a unique solution to the system. Although the majority of the parameters correspond to a well-defined minimum of the cost function, a few are associated with a plateau at the minimum error value (see Supplementary text, Parameter sensitivity, for details). This is not surprising or uncommon behavior in biological dynamical systems, as they may often be decomposed into fast and slow subsystems. In such cases, when a parameter that governs kinetics of the fast subsystem is made sufficiently large, a wide range of parameter values will yield similar dynamical behavior. This behavior indicates that there is a threshold value at which the kinetics of the fast subsystem controlled by this parameter effectively becomes instantaneous from the temporal perspective of the slower subsystems, which are rate limiting. Such robustness to parameter variation is not uncommon in biological systems. Thus, the final parameter values of the model presented here are not necessarily intended to quantitatively capture and reproduce the kinetic features of all the underlying biophysical processes associated with each gating transition of the channel protein, but rather as appropriate values representing a simplified mechanistic model system capable of reproducing experimentally observed channel behavior. Finally, a few commonly used LTCC models do not assume independence between channel activation and CDI (Greenstein and Winslow, 2002; Faber et al., 2007; Miriyala et al., 2008; Decker et al., 2009). Thus, integrating our CDI formulations into these LTCC models may be nontrivial in that the combined models would require a large number of states (∼80 with a five-state activation model). However, the model size is not a limitation for stochastic gating simulations, and methods for handling such models in ODE-based simulations have been established previously (Greenstein et al., 2006).

Overall, beyond predicting the CDI profile of LQTS-associated calmodulinopathies, our model could be applied in a broader context to gain insights into the (patho)physiology of cardiac electrical activity. Previous studies have implied that CDI, compared with VDI, may be the dominant form of feedback regulation of Ca_V1.2 in the context of the cardiac action potential (Sun et al., 1997; Morotti et al., 2012) and the development of arrhythmia (Tanskanen et al., 2005; Morotti et al., 2012). Therefore, it is critical to use a realistic mathematical representation of CDI that can capture the distinct kinetics and cooperativity of the two lobes of CaM. Built from experimentally derived knowledge of fundamental biophysical gating mechanisms, this bilobal kinetic model of CDI can be used in cell and tissue scale models to make reliable predictions about cardiac electrical activity. Finally, the model can be used within other cellular contexts, including neurons, vascular and gastrointestinal smooth muscles, and chromaffin cells (Zamponi et al., 2015). Therefore, this bilobal CDI model could be an invaluable asset to the field of computational electrophysiology.

David T. Yue, who originally conceptualized this work, passed away on December 23, 2014. His mentorship, wisdom, and kindness are greatly missed. We thank Dr. Shin Rong Lee for advice and training with the photolysis system and thank Wanjun Yang for dedicated technical support. We also thank Dr. Gordon Tomaselli for providing valuable advice and discussions, and Dr. Manu Ben-Johny and members of the Calcium Signals Lab for ongoing feedback and support.

This work was supported by grants: American Heart Association Pre-doctoral Fellowship (W.B. Limpitikul) and National Institutes of Health R01HL105239 (J.L. Greenstein and R.L. Winslow) and R01MH065531 (W.B. Limpitikul, D.T. Yue, and I.E. Dick).

The authors declare no competing financial interests.

Author contributions: W.B. Limpitikul collected and analyzed data and wrote and edited the manuscript. J.L. Greenstein, R.L. Winslow, and I.E. Dick supervised the overall project and wrote and edited the manucript. D.T. Yue originally envisioned and conceptualized this work.

Richard W. Aldrich served as editor.