Signaling proteins such as ion channels largely exist in two functional forms, corresponding to the active and resting states, connected by multiple intermediates. Multiparametric kinetic models based on sophisticated electrophysiological experiments have been devised to identify molecular interactions of these conformational transitions. However, this approach is arduous and is not suitable for large-scale perturbation analysis of interaction pathways. Recently, we described a model-free method to obtain the net free energy of activation in voltage- and ligand-activated ion channels. Here we extend this approach to estimate pairwise interaction energies of side chains that contribute to gating transitions. Our approach, which we call generalized interaction-energy analysis (GIA), combines median voltage estimates obtained from charge-voltage curves with mutant cycle analysis to ascertain the strengths of pairwise interactions. We show that, for a system with an arbitrary gating scheme, the nonadditive contributions of amino acid pairs to the net free energy of activation can be computed in a self-consistent manner. Numerical analyses of sequential and allosteric models of channel activation also show that this approach can measure energetic nonadditivities even when perturbations affect multiple transitions. To demonstrate the experimental application of this method, we reevaluated the interaction energies of six previously described long-range interactors in the Shaker potassium channel. Our approach offers the ability to generate detailed interaction energy maps in voltage- and ligand-activated ion channels and can be extended to any force-driven system as long as associated “displacement” can be measured.
Allosteric control of enzymatic activity is a well-known phenomenon in biological systems (Goodey and Benkovic, 2008). It often involves a ligand or metabolite that has no stereochemical resemblance to either products or substrates but exerts control on protein function by binding to a distal regulatory site. However, allosteric regulation is not just limited to substances that exert chemical potential. Activity of proteins such as ion channels can be controlled by voltage, mechanical stretch, or even heat (Hille, 2001). Understanding the molecular mechanisms underlying energy transduction in proteins requires identification of specific pathways involved in transmission of information from sensory modules to catalytic centers. Structural methods like nuclear magnetic resonance and x-ray crystallography provide high-resolution snapshots of conformational changes that underpin molecular communication, but ultimately direct measurement of energetic coupling is also necessary.
Interaction energies between specific residues involved in protein folding can be determined by way of mutant cycle analyses (Fig. 1 A; Ackers and Smith, 1985; Horovitz and Fersht, 1990; Horovitz et al., 1994; Di Cera, 1998). In this approach, perturbation energies are evaluated when a specific site (say X) is perturbed in the native protein and in the background of a secondary perturbation (say Y). If the two perturbation energies are equal, then it implies that the two sites are energetically independent, whereas unequal perturbation energies imply that the two sites interact. This methodology was first used to study protein–protein interactions in tRNA synthetase (Carter et al., 1984) and has now seen widespread application in the study of protein folding and conformational changes (Serrano et al., 1990; Schreiber and Fersht, 1995; Ranganathan et al., 1996). Of central importance is this question: How is the perturbation energy in each step evaluated?
In the case of a binding reaction such as those involving protein–protein or protein–ligand interaction, free-energy changes associated with perturbations can be evaluated by directly monitoring binding. In many cases, however, functional activity is used as a surrogate measure of free energy of conformational transition. For instance, free energy of activation of voltage-gated and ligand-activated ion channels is evaluated based on their functional responses as evidenced in conductance-voltage (G-V) or conductance-ligand (typically referred to as dose response) relationships. If the channel has a single conducting state, these relationships reflect the channel open probabilities (PO) at different voltages or ligand concentrations and are fitted to a logistic function. For voltage-activated processes, the free energy of activation of the channel is estimated as ΔGapp = zFV1/2 (V1/2 is the voltage that elicits half-maximal response, and z reflects the number of charges transferred during channel activation; Fig. 1 B). These free-energy measures can be combined with double mutant cycle analyses to estimate pairwise interaction energies between various sites in a voltage (Yifrach and MacKinnon, 2002; Sadovsky and Yifrach, 2007; DeCaen et al., 2008, 2009, 2011; Zandany et al., 2008; Yifrach et al., 2009; Wall-Lacelle et al., 2011; Cheng et al., 2013; Chamberlin et al., 2014)- or ligand-activated ion channels (Kash et al., 2003; Gleitsman et al., 2009; Shanata et al., 2012). This approach will be hereby referred to as the functional mutant cycle (FMC).
In most known voltage-gated ion channels, however, the process of channel activation comprises multiple intermediate closed states, before the channels finally open (Bezanilla et al., 1994; Zagotta et al., 1994a; Schoppa and Sigworth, 1998). The occupancies of the ensemble of nonconducting states are not adequately captured by the G-V curves, and their energetic contributions remain indeterminate. More significantly, the uncertainty associated with interaction energies quantified using simplified or empirical free-energy measures prevents a direct comparison with molecular simulations based on high-resolution structures. In some instances, multistate kinetic models are fitted to experimental behavior of proteins and the model parameters are directly used for mutant cycle analyses (Lee and Sine, 2005; Gupta and Auerbach, 2011). Although such an approach is likely to be better than simply measuring G-V curves, its use is limited to systems whose gating behavior can be adequately described by a well-constrained gating scheme. We should note that dose–response curves and G-V curves in allosteric systems under limiting conditions can lend themselves to linkage analysis to probe interaction energies (Chowdhury and Chanda, 2010; Sigg, 2013). However, this approach is also extremely time consuming and cannot be easily extended to obligatorily coupled systems (Chowdhury and Chanda, 2012b).
Recently, we described an alternate methodology to estimate the voltage- and ligand-dependent change in free energy during channel activation (Chowdhury and Chanda, 2012a, 2013), which was an extension of linkage analysis developed by Wyman and Gill (1990). We showed that this net free-energy change (ΔGnet) for a voltage-dependent channel can be calculated by measuring median voltage of activation from gating charge voltage curves (Fig. 1 C). An advantage of this approach is that it does not require us to build multistate kinetic models to calculate free energy of activation (Miller, 2012). Here, we test the proposition that these free-energy measurements can be combined with mutant cycle analyses to determine interaction energies between sites in a model-independent fashion. Previous studies on the Shaker KV channel identified several putative interactors in the pore domain that were proposed to impact the last concerted pore opening transition (Yifrach and MacKinnon, 2002; Sadovsky and Yifrach, 2007). We find that the majority of these interaction pairs do not contribute significantly to the overall free energy of activation of the channel. Only a single residue pair, A391-T469, exhibits significant energetic nonadditivity, which in all likelihood reflects a long-range interaction between them mediated by networks of intervening residues. Extensive numerical simulations, using different multistate gating models, clearly show that generalized interaction-energy analysis (GIA) yields a self-consistent estimate of energetic nonadditivities. We expect that this approach will allow us to pursue large-scale analysis of interaction networks that underlie gating transitions in voltage- and ligand-gated ion channels.
MATERIALS AND METHODS
Mutagenesis and expression in Xenopus laevis oocytes
All clones used in this study were derived from a cDNA of the inactivation-removed Shaker KV channel (Δ6–46) cloned into the pBSTA vector. Mutations were introduced by PCR using mismatch mutagenic primers (QuikChange; Agilent Technologies). All mutations were confirmed by sequencing the whole cDNA. For gating current measurements, mutations were introduced into the background of the W434F mutation and rendered the channel nonconducting without compromising other aspects of gating (Perozo et al., 1993). Mutant cDNAs were linearized using a NotI enzyme (New England Biolabs, Inc.) and transcribed into cRNAs using mMESSAGE mMACHINE T7 kit (Life Technologies).
Xenopus oocytes were removed surgically and treated with 1 mg/ml collagenase for 1–1.5 h to remove the follicular layer. Before injection, oocytes were incubated in ND-96 solution supplemented with 100 µg/ml gentamicin at 18°C, typically for 12–24 h. 50 nl cRNA at a concentration of 50–200 ng/µl was injected into oocytes. After injection, the oocytes were kept at 18°C in a solution containing 100 mM NaCl, 2 mM KCl, 1.8 mM CaCl2, 1 mM MgCl2, 5 mM HEPES, 0.1 mM DTT, and 0.2 mM EDTA, supplemented with 100 µg/ml gentamicin and 100 mg/ml bovine serum albumin. Ionic current measurements were performed 12–24 h after injection, whereas gating currents were measured 2–7 d after injection.
Ionic currents were measured on a cut-open oocyte voltage clamp (COVC) set-up (CA-1B; Dagan Corporation) as described previously (Gagnon and Bezanilla, 2010). The external solution used was 105 mM NMG-MES (N-methyl-d-glucamine methanesulfonate), 10 mM K-MES, 2 mM Ca-MES, and 10 mM HEPES, pH 7.4. The internal solution was 115 mM K-MES, 2 mM EGTA, and 10 mM HEPES, pH 7.4. Gating currents were measured either on a COVC or two-electrode voltage (TEV) clamp set-up. Some of the mutants investigated in this study exhibited low expression and their gating currents were measured on a TEV clamp set-up. The external solution used for gating current measurements in both set-ups was 115 mM NMG-MES, 2 mM Ca-MES, and 10 mM HEPES, pH 7.4. The internal solution used for gating current measurements on the COVC set-up was 115 mM NMG-MES, 2 mM EGTA, and 10 mM HEPES, pH 7.4. The recording pipette resistance for all electrophysiological measurements was 0.2–0.5 MΩ. Analogue signals were sampled at 20–250 kHz with a Digidata 1440 or 1320 interface (Molecular Devices) and low-pass filtered at 10 kHz.
Ionic currents were obtained by applying 100-ms-long depolarizing pulses from −120 to 80 mV in 2.5- or 5-mV increments. The holding potential was −120 mV. Capacitive transients and linear leak currents were subtracted online using the P/−4 method, during which the holding potential was −120 mV. After baseline subtraction, peak tail current amplitudes, elicited by repolarization pulses to −120 mV, were used to generate the relative open probability versus voltage (I/Imax) curves. Gating currents were obtained by applying a 50-ms-long depolarizing pulse to voltages from −120 to 20 mV (in 5-mV intervals). For measurements using COVC, the holding potential used was −120 mV, whereas on the TEV, the holding potential was −90 mV as it was not possible to hold the oocytes at −120 mV. Depolarization pulses were preceded and followed by 50-ms pre- and postpulses to −120 mV. Q-V curves for the W434F mutant measured in either of the two setups were superimposable. The capacitive transient and linear leak currents were subtracted online using the P/−4 or P/−8 method, with a subsweep holding potential of −120 or −90 mV (on COVC or TEV, respectively). After baseline readjustments, the on-gating current records were integrated over the duration of the depolarization pulse to obtain the gating charge displaced, which was used to compute the fractional gating charge displacement versus V curve (Q/Qmax vs. V or Q-V).
The I/Imax curve for each mutant was obtained by averaging the curves obtained from three to six oocytes. The curve was fitted to the Boltzmann equation:
where zapp is the Boltzmann slope and V1/2 is the voltage that elicits half-maximal response. The Boltzmann measure of free-energy change, ΔGapp, for each mutant was calculated as ΔGapp = zappFV1/2. The uncertainty associated with ΔGapp estimation (δΔGapp) was calculated as
where δzapp and δV1/2 are the standard error associated with estimation of zapp and V1/2, respectively.
The nonadditivity in an FMC, ΔΔGFMC, was calculated as
where the subscripts indicate the WT channel, double mutant (S12), or the two single mutants (S1 or S2). The uncertainty associated with ΔΔGFMC (δΔΔGFMC) was calculated as
where the terms in the brackets indicate the uncertainty associated with the ΔGapp estimates of the WT and double and two single mutants.
The fractional gating charge displacement curves for all of the mutants were obtained by averaging measurements performed on three to six oocytes. The median voltage of activation, VM, for each Q-V curve was extracted by calculating the area between the Q-V curve and the ordinate axis, using the trapezoid method. For a Q-V curve with n points, the VM is calculated as
where Qi and Vi are the ith point on the QfV curve. The net free energy of activation of the channel is calculated as ΔGC = QmaxFVM, where Qmax is the maximum number of charges transferred during voltage-dependent activation of the channel. For all of our calculations, we used a Qmax of 13.2 (Schoppa et al., 1992; Aggarwal and MacKinnon, 1996; Seoh et al., 1996). Although, Qmax for each of the mutations was not measured individually, it is unlikely that any of the mutants studied in this paper alter Qmax as they are not the primary gating charge–determining residues of the channel (Aggarwal and MacKinnon, 1996; Seoh et al., 1996). The uncertainty in ΔGC was calculated as QmaxFδVM, where δVM is the standard error of the VM estimation.
The nonadditivity in a mutant cycle analysis was calculated using the median measure of free-energy change; this nonadditivity, ΔΔGGIA, was calculated as
The standard error associated with ΔΔGGIA (δΔΔGGIA) was calculated as
where δ(VM)WT, δ(VM)S1, δ(VM)S2, and δ(VM)S12 are the uncertainties (standard error of the mean) associated with VM measurement of the WT channel and the single and double mutant channels, respectively.
All simulations were performed with MATLAB version R2012b. The nonadditive energies for randomly generated mutant cycles were performed as follows. First we assume that all of the three mutants and the WT channel constituting the mutant cycle undergo voltage-dependent activation following the same scheme. This assumption implies that there is at least one discrete state Markov model (gating scheme) that can adequately describe the gating properties of all the mutants and the WT channel, constituting the mutant cycle. Various constructs of the cycle differ in the values of model parameters, although the number of closed and open states and their connectivities remain the same. This assumption is invoked only to simplify the algebraic representation of ΔGnet and ΔΔGtrue (in terms of the different equilibrium parameters of the gating scheme) when we compare energetic nonadditivities using GIA and FMC. Each gating scheme is described by n equilibrium parameters whose magnitude at 0 mV are K1, K2,…, Kn and all of which have a standard exponential voltage dependence determined by parameters z1, z2,…, zn, respectively. For the ith equilibrium constant, the range of values (i.e., the parameter space) is determined to be to and similarly the voltage-dependent parameters also have a defined range of values to In each mutant cycle, the model parameters for the WT were fixed, but those of the mutants were selected randomly within the given parameter space. For each mutant in the cycle, “2n” random numbers between 0 and 1 (a uniform random number generating function was used) were generated, each of which specified the extent of perturbation on the equilibrium constants or the voltage dependence parameters. The equilibrium constants of the mutant was assigned as where R is the random number, whereas the voltage-dependent parameters were assigned as (the random numbers used for perturbation of Ki and zi were independently generated). This process was repeated for all three mutants of the cycle and represents a situation where single and double mutants have multiple effects on channel energetics. The PO-V and Q-V curves of the mutants were generated from which ΔGapp and ΔGnet was extracted. These parameters were then used to compute ΔΔGFMC and ΔΔGGIA for the cycle. This process was repeated for >600 cycles in two instances: first where the channel and its mutants follow the ZHA activation scheme (see Fig. 5 A) and second where they follow MWC activation scheme (see Fig. 6 C). This strategy allowed us to test the accuracy of ΔΔGGIA over a large parameter space.
The parameter range sampled for the simulations performed using the ZHA scheme was [0.005–5,000], [0.1–100,000], L0: [0.01–10000], z1: [0.5–3.5], z2: [0.5–3.5], and zL: [0.25–2.5]. The parameters for the WT channel were fixed at 5, 100, L0: 10, z1: 2, z2: 1.5, and zL: 0.5 (which are close to the parameters for the Shaker KV channel reported previously [Zagotta et al., 1994a; Ledwell and Aldrich, 1999]). The parameter space sampled for the simulations performed using the MWC scheme was J0: [0.00001–10], L0: [10−9–10−3], D: [1–100], zJ: [0.1–4], and zL: [0.1–2] (D was maintained voltage independent in all cases). The parameters for the WT channel were fixed at J0: 0.03, L0: 10−6, D: 25, zJ: 0.6, and zL: 0.3 (which are close to the parameters for the BK channel model reported previously [Horrigan and Aldrich, 2002]).
Online supplemental material
Figs. S1 and S2 show a family of gating and ionic current traces for Shaker potassium channel mutants that were investigated in this study.
Principle of the GIA
Let us consider a protein that exists in a passive form, S1, which upon action of an external stimulus (such as ligand, voltage, etc.) undergoes a series of conformational changes to reach its final active form, Sn+1, via n sequential transitions (involving n − 1 intermediate states, Si, i = 2, 3,…, n − 1). The free energy of each conformational state is Gi. On such a system we implement the mutant cycle analysis (Figs. 1 A and 2) to deduce the interaction energy between two residues, X and Y.
The residue X is perturbed in the native protein and in the background of the secondary perturbation, Y. As a consequence of these perturbations, the free energies of multiple intermediates may change. The interaction energy between the two residues in the state Si (i = 1, 2,…, n), can be described by the equation
where the superscripts on the free-energy terms on the right side of the equation indicate the unperturbed and singly and doubly perturbed systems, as shown in Fig. 2. Equations, analogous to Eq. 1, may be written for each of the conformational states of the protein.
To understand how these interactions drive gating transitions, we would like to measure the changes in these interaction energies when the channel undergoes a conformational change. From an experimental standpoint, we can measure the free-energy change associated with a conformational change, say between states Si and Si+1, by measuring the equilibrium constant, Ki, for the transition. This allows us to evaluate the difference in the interaction energies between the two states in two different conformations as
which may be rewritten as
In Eq. 2, ΔΔGi reflects and ΔGi→i+1 is the free-energy difference between the states Si and Si+1, which can be written as −RTlnKi. Ωi reflects the nonadditivity of the equilibrium constants in the mutant cycle.
Next, suppose the nonadditive perturbation energy associated with each of the n conformational transitions (ΔΔGi, i = 1, 2,…, n) are evaluated and we sum all of the ΔΔGi measures to obtain the net nonadditivity of the two perturbations, ΔΔGnet:
where ΔGi→i+1 is the free-energy difference between the states S1 and Sn+1. This net energetic nonadditivity reflects the difference in the interaction energies between X and Y in the initial passive conformation and final active conformation. A protein might transit between its two limiting states following multiple pathways. Because most biological macromolecules obey the principle of microscopic reversibility, the net free-energy changes across all such pathways are necessarily identical. This implies that ΔΔGnet is path independent, and thus Eq. 3 holds for systems undergoing multistate transitions in modes that are sequential or parallel or combinations thereof.
How can we extract ΔΔGnet from experimental measurements? Recently, we described an approach to extract the ΔGnet for conformational change in proteins, driven by an external stimulus, in a model-independent manner (Chowdhury and Chanda, 2012a, 2013). The approach involves measuring the “conjugate displacement” associated with the stimulus, followed by an integral transformation of the conjugate displacement versus “stimulus intensity,” which directly yields the total work done during conformational change in the protein (Chowdhury and Chanda, 2012a, 2013; Sigg, 2013). This method was described in detail for voltage-gated ion channels, for which the conjugate displacement curve is the gating charge displacement versus voltage (Q-V) curve. From the measured Q-V curves we extract the median voltage of activation, VM, and evaluate ΔGnet as QmaxFVM, where Qmax is the maximum amount of gating charges transferred during channel activation (or in other words the charge per channel). By incorporating such a free-energy measure in Eq. 3, we can obtain a measure of ΔΔGnet (Eq. 3). Because this approach of combining the mutant cycles with measurements of conjugate displacement curves is referred to as GIA, the ΔΔGnet in Eq. 3 will be renamed ΔΔGGIA.
Experimental comparison of FMC analysis and GIA
FMC analysis and GIA use orthogonal experimental measurements to estimate the interaction energies between two sites. How do the interaction energies computed by these two methods compare against each other? In the prototypical Shaker KV channel, previous FMC analysis has identified several strongly interacting residues (Yifrach and MacKinnon, 2002; Sadovsky and Yifrach, 2007). Some of the interaction partners are >10 Å away in the protein structure, which has led to the view that voltage-dependent channel opening involves dynamic rearrangements of long-range interaction networks. To estimate these interaction energies accurately, we reevaluated them using the GIA approach.
We focused on four residues (A391, E395, T469, and V476) located in the pore domain of the channel, a region likely to undergo significant conformational changes during channel gating (Fig. 3 A; Yifrach and MacKinnon, 2002; Sadovsky and Yifrach, 2007). First, we applied FMC to measure the interaction energy between sites E395 and T469. Two single and one double mutant were generated by mutating the sites to Ala, and their relative PO-V curves were measured (along with that of the WT channel) by tail current analysis (Fig. 3 B and Fig. S1). Both single mutations result in large shifts in the PO-V curves, and consistent with a previous study (Yifrach and MacKinnon, 2002), FMC analysis shows a large nonadditivity (ΔΔGFMC ∼7.5 kcal) of the two perturbations (Table 1). Next, for each of the constructs, Q-V curves (Fig. 3 C and Fig. S1) were measured, which allowed us to calculate perturbation energies caused by mutations by using median analysis (Table 2; Chowdhury and Chanda, 2012a). Strikingly, GIA showed that the perturbations at these sites are energetically independent (ΔΔGGIA ≈ 0; Table 3).
Next, we extended this comparison to other residue pairs that have been previously proposed to be involved in interactions (Fig. 4, Fig. S2, and Tables 2 and 3). Q-V curves of single mutants (A391V, E395A, T469A, and V476A) and all pairwise combination (double) mutants were measured to calculate ΔΔGGIA. Studies in proteins show that interaction energies between noncharged residues range between 0.5 and 1.0 kcal/mol (Horovitz, 1996). We set the cut-off for interaction in a single subunit at 0.45 kcal/mol, and therefore, for a tetrameric channel it will be 1.8 kcal/mol. This is also above our experimental error associated with ΔΔGGIA estimates, which was ∼0.6 kcal (≈RT). In all of the six possible pairs, we found that that ΔΔGGIA and ΔΔGFMC (reported previously) do not agree with each other. More importantly, except for the AT pair (i.e., A391-T469), all of the pairwise nonadditivities evaluated by GIA were below the significance level, suggesting that they are energetically independent.
Numerical simulations of GIA and FMC
This discrepancy in interaction energies evaluated using FMC and GIA prompted us to examine the robustness of the two approaches in measuring energetic nonadditivities. We performed numerical simulations using the 16-state ZHA model (Fig. 5; Hoshi et al., 1994; Zagotta et al., 1994a,b; Ledwell and Aldrich, 1999), which describes the voltage-dependent activation of the Shaker KV channel. In this model, voltage sensor activation is described as a two-step sequential process occurring independently in different subunits, and once all voltage sensors are activated, a final concerted transition opens the channel pore. We constructed a hypothetical mutant cycle in which two mutations were envisioned to affect only the equilibrium constant of the last concerted transition and the double mutant to affect the same transition, additively (Fig. 5 A). Thus, by design, the cycle ensures that the two perturbations are energetically independent, and hence nonadditivity measurements should yield a null result. While holding the other equilibrium parameters constant, the magnitude of perturbation of each mutant was varied over 12 orders of magnitude, and in each case we calculated the ΔΔG using FMC (by simulating the PO-V curves) and GIA (by simulating the Q-V curves). As shown in Fig. 5 (B and C), the FMC simulations show strong nonadditive energies, whereas the GIA simulations show that the nonadditivity is zero in all cases. This simulation suggests that even under limiting conditions when perturbations affect a single transition of the gating scheme, FMC analysis may report nonreal interaction energies.
Next, we considered a more general scenario where single mutations affect multiple transitions randomly and the double mutants may affect these transitions either additively or nonadditively. A hypothetical mutant cycle was created wherein the WT reference channel and the three mutants constituting the cycle were envisioned to be gated via the ZHA model. For each of the three mutants comprising a thermodynamic cycle, the values of the equilibrium parameters were randomly chosen from within a parameter space (see Materials and methods). The PO-V and Q-V curves in each case were simulated and the ΔΔGFMC and ΔΔGGIA for the cycle were computed. This process was repeated for a large number (>600) of mutant cycles, and ΔΔGFMC and ΔΔGGIA for each cycle were then compared against the true nonadditive energy ΔΔGtrue (Fig. 6, A and B), derived directly from the model equilibrium constants. ΔΔGtrue was calculated as
where the superscripts WT, S12, S1, and S2 indicate the equilibrium constants for the WT channel and the double and the two single mutant channels, respectively. The comparison reveals that ΔΔGGIA and ΔΔGtrue are identical for all cases (Fig. 6 B), whereas ΔΔGFMC is not correlated with ΔΔGtrue (Fig. 6 A).
Similar simulations were performed for an MWC-type allosteric scheme (Fig. 6 C). According to this scheme, voltage sensor activation and pore opening represent preexisting equilibria wherein the pore is intrinsically biased toward the closed state and voltage sensor activation causes a shift in this bias toward the open state through allosteric interactions. Such an activation scheme has been shown to accurately describe the voltage-dependent gating of the BK channel (Horrigan and Aldrich, 1999, 2002). For >600 mutant cycles generated by a random sampling strategy, we obtained ΔΔGFMC and ΔΔGGIA, whereas ΔΔGtrue was calculated as
where J and L represent the intrinsic equilibrium constant of activation, at 0 mV, for the voltage sensor and the pore and D represents the allosteric interaction factor. We find that ΔΔGFMC and ΔΔGtrue are highly divergent (Fig. 6 D), whereas ΔΔGGIA and ΔΔGtrue are equal in all tested cases (Fig. 6 E). The discrepancy between ΔΔGGIA and ΔΔGtrue is thus not just limited to the systems activating via the ZHA model but also extends to other multistate voltage-dependent systems (see Appendix for additional examples). These findings establish that GIA but not FMC can report energetic nonadditivities (residue-specific interaction energies) in a self-consistent manner irrespective of whether the perturbations affect single or multiple transitions.
In this study, we describe an experimental methodology to quantitatively determine the contribution of interaction energies between two sites of a protein to the overall free-energy change associated with a stimulus-driven conformational change of the protein. Our method is based on determining the energetic consequence of a perturbation in the presence of a secondary perturbation, by measuring the conjugate displacement of the process for each of the perturbed systems. This is fundamentally different from the canonical mutant cycle analyses, wherein functional activity of a protein is measured and used to empirically quantify the free energies. We show that the free energies of perturbations calculated via a median transformation of the conjugate displacement curve allow us to directly calculate the interaction energies between specific sites of a multistate protein, in a self-consistent manner, even in instances when the perturbations affect multiple transitions of the protein.
We implemented GIA on the voltage-gated Shaker KV channel, for which the conjugate displacements associated with the voltage-dependent transformation can be conveniently obtained by measuring the gating charge displacement versus voltage (Q-V) curves. We used GIA to calculate the interaction energies between several pairs of residues on the pore domain of the channel, which undergoes large conformational changes during channel gating (Yellen, 1998). Many of the perturbed pair of sites are several angstroms away in the structure, and previous FMC studies have suggested that they constitute a long-range energetically coupled network of residues that are crucial for channel gating (Yifrach and MacKinnon, 2002; Sadovsky and Yifrach, 2007). However, using GIA, we find that all except for one (A391-T469) were energetically independent of each other. We should note that some of these differences in interaction energy measurements could arise because of the use of different backgrounds. We use W434F nonconducting mutant in our study, whereas conducting Shaker potassium channel was used as a background for FMC studies (to be discussed later). Numerical simulations of allosteric (MWC type) and the quasi-sequential nonallosteric (ZHA type) models show that GIA can provide estimates of nonadditivities over a large parameter space even when mutations affect different transitions during a multistep activation process.
In this study, the GIA approach is implemented by introducing alanine substitutions at each of the test sites, which removes all atoms in the native side-chains except the Cβ atom. Thus, the deduced energetic nonadditivity reflects the excess contribution of native residue pairs to the overall channel energetics, relative to the double alanine pair at the respective sites. Although a glycine substitution might be deemed better than alanine, the backbone flexibility introduced by glycine substitutions could introduce structural distortions to the protein and thus is not preferable in most instances (Faiman and Horovitz, 1996; Di Cera, 1998). Additionally, native alanine residues are frequently mutated to valine, which should not be directly compared with alanine-based mutant cycles (Yifrach and MacKinnon, 2002).
Interpreting GIA and FMC in terms of energy landscape
The false positives that are observed in the FMC approach can be rationalized by considering a simple three-state model of a channel with two closed states (C0 and C1) and one open state (O; Fig. 7). A G-V curve samples the occupancy of only the O state. As a consequence, the G-V–based free-energy estimate is dominated by the energy difference between the open state and most stable closed state. Mutations that alter the stability of the intermediate closed state or states could change the reference closed state. In the example shown, G-V–based energetics are informed by the energy difference between C1 and O for the WT but C0 and O for the other three mutants. Because of these differences in reference states, G-V–based perturbation energies are not comparable across mutations, which may contribute to inaccurate interaction energy estimates (Fig. 7). This problem is exacerbated as the number of intermediate closed states increase. In contrast, GIA always measures the energy difference between first and last states. This ensures that in the given example the Q-V–based energies always reflect the energy difference between the states C0 and O. In other words, the GIA approach allows us to reduce the complexity of the multistate gating process to a simple comparison between binary states.
We note that energetic additivity in GIA will arise in three instances: (1) the sites truly do not interact, (2) the sites interact only in the intermediate states, or (3) the interaction between the two sites does not change between the initial resting and final activated states. Thus, additivity is not a proof for absence of interaction, but when nonadditivity is observed, one can be certain that the residue pairs synergistically contribute to the gating process.
When channel activation is truly a two-state process, the nonadditive energies of FMC will be self-consistent (Horovitz and Fersht, 1990). Another scenario where FMC might be self-consistent is when the Q-V and relative PO-V curves of the channel and the mutants are widely separated along the voltage axis. In such instances, the PO-V curve can be approximated to describe a single-step transition between the last closed and open states, which is reminiscent of specific mutants of the Shaker KV channel (Ledwell and Aldrich, 1999; Soler-Llavina et al., 2006). ΔΔGFMC obtained from Boltzmann fits to the PO-V curves would reflect the nonadditivity in the last step of channel gating; however, the nonadditive energies in the earlier steps or the overall gating process will not be accessible. In contrast, as shown earlier (Eq. 3), ΔΔGGIA reflects the summation of energetic nonadditivities in all of the steps from the initial resting state of the channel to the final open state of the channel. Irrespective of whether the nonadditivity arises because of changes in interresidue interactions in the initial, intermediate, or final steps, ΔΔGGIA will reflect this nonadditivity as long as energetic nonadditivity in one transition is not compensated by that in another transition. Identifying the specific transitions that contribute to nonzero ΔΔGGIA will require additional kinetic analysis.
Considerations for gating charge measurements
In the Shaker KV channels, gating charge displacement curves have been measured in many different ways. In this study, we used the nonconducting W434F mutant to measure the Q-V curves. This mutant eliminates ion conduction by accelerating C-type inactivation, resulting in a putatively collapsed outer pore (Yang et al., 1997; Cordero-Morales et al., 2011). Gating charge of different channel mutants can also be obtained in the background of a second nonconducting mutant of the Shaker KV channels: V478W, which stabilizes a hydrophobic seal responsible for pore closure and thereby occludes ion flux (Kitaguchi et al., 2004). A third strategy to measure gating currents is to use a high-affinity peptide toxin, such as Agitoxin, which blocks ion conduction (Aggarwal and MacKinnon, 1996). The final alternative is to measure gating currents under conditions when all permeant ions are washed off. Although the Q-V curves obtained from each of the four strategies are generally considered to be equivalent, there may be subtle differences between them, especially for mutant channels. Therefore, our interaction energies measured here should be considered as those obtained in the background of W434F mutant. A thorough and careful study will be needed to assess and compare these different strategies to extract the true perturbation energies for the WT channel.
Another consideration is whether VM or V1/2 of the Q-V curves can be used interchangeably for free-energy calculations. Although the Q-V curves of Shaker KV channel and its mutants exhibit a slight asymmetry, in most instances, the Q-V curves can be adequately fitted to the symmetric Boltzmann equation and the median voltage of charge transfer (VM) will be often similar to V1/2, the voltage at which half of the gating charges have moved. In such situations, Boltzmann fit–derived V1/2 values can be used instead of VM values to compute the energetic nonadditivities despite the fact that the underlying gating is seldom restricted to two states. However, when the Q-V curves are split, V1/2 values depend on the choice of model used for fitting and thus might not be unique (Chowdhury and Chanda, 2012a). In contrast, VM is not obtained through model-specific equations and is characteristic for each Q-V curve.
Implementation of GIA requires prior knowledge of Qmax, the maximum gating charge transferred during full-scale activation, of the channel, and its mutants. Measurement of Qmax for a channel is a nontrivial problem and cannot be estimated simply from the slope of the Q-V curve (Bezanilla and Villalba-Galea, 2013). However, for the system under investigation, the gating charge–determining residues reside in the voltage sensor and are highly specific (Ahern and Horn, 2004). The specific sites that have been investigated in this study are uncharged, lie outside the electric field, and when perturbed are unlikely to significantly alter the Qmax of the channel. However, in implementing GIA to test for interactions between residues in the voltage sensor, it might be important to calibrate Qmax of the mutant channels to establish that the Qmax is not altered significantly by the perturbations.
Considerations for gating involving voltage-independent transitions
For many ion channels the maximum open probability () at depolarized potentials does not reach unity because the transition between the last closed and open states is voltage independent. In such situations, a correction factor () needs to be added to QmaxFVM to obtain an accurate measure of the free energy of channel activation (Chowdhury and Chanda, 2012a). Incorporating this correction factor, the energetic nonadditivity in GIA is
In Eq. 4, the second term on the right side of the equation incorporates the maximum open probabilities of the WT and single (S1 and S2) and double (S12) mutants, as indicated in the superscript. The magnitude of this correction factor becomes significant (>1.8 kcal) only when becomes larger than ∼25, which would imply that of at least one of the mutants is 25-fold lower than that of the WT channel. Although such an effect is possible, for many Shaker KV channel mutants (in the voltage sensor and the pore), has been measured and found to be similar to that of the WT channel (Seoh et al., 1996; Ding and Horn, 2002). This implies that although measurements will improve the accuracy of the free-energy estimate, its contribution is likely to be small enough so as to not interfere with the identification of strongly interacting residue pairs. A similar consideration would apply to channels that exhibit constitutive opening even under hyperpolarizing conditions. As noted previously (Chowdhury and Chanda, 2012a, 2013), under situations when the of the channel is significantly large, a correction factor amounting to needs to be incorporated in the overall free-energy equation to accurately quantify the free-energy difference between the initial and final state (along the charge coordinate) of the channel. However, such a correction factor will be significant (at level of the thermal energy) only when > 0.7, which, although possible, is likely to be a rare occurrence.
We should emphasize that GIA provides the net free energy of interresidue interaction associated with gating transitions that are driven by externally applied force. The corrections described above are to take into account the free-energy contributions of transitions that are not driven by external force but are still functionally relevant. This means that GIA without any corrections provides us a model-free method for obtaining interaction energy change going from the initial force-dependent state to the final force-dependent state.
Finally, as a note of caution, we should add that the interaction energies measured by GIA are true interaction energies only if the free energy of the system can be simply expressed as a sum of components stemming from specific interactions (Mark and van Gunsteren, 1994). This underlying assumption may not hold true especially when higher-order mutant cycles are constructed, and therefore these interaction energies should be regarded as empirical interaction energies. But this assumption underpins all thermodynamic cycle analyses, and unlike FMC measurements, GIA measurements of interaction energies are self-consistent with multistate gating schemes.
Application of FMC analyses to multistate systems such as ion channels have long relied on empirical free-energy metrics that have been argued to be correlated (occasionally linearly) with the true free-energy change in the system. Such empiricism, although occasionally necessary, is associated with significant uncertainty and, as demonstrated through the simulations performed in this study, can potentially obscure essential molecular and thermodynamic features of the system. The firm thermodynamic foundation of GIA greatly reduces such ambiguities and therefore offers an indubitable strategy to identify residue-level interactions that account for the thermodynamic changes associated with structural dynamics of a protein as it is driven by an external force. Although the application of GIA based on normalized conjugated displacement curves needs to be carefully implemented, as has been described in detail in the previous sections, it offers a clear advancement over the traditional FMC analyses by circumventing the necessity of a two-state process approximation.
Although our discussion in this manuscript is focused mainly on the application of GIA for voltage-dependent systems, this principle can be extended to other stimulus-driven systems. For a ligand-driven transformation, the direct ligand-binding curve is the conjugate displacement curve, and by measuring the median ligand activity, we will be able to estimate the total free-energy change associated with the ligand-driven transformation in a model-independent fashion (Chowdhury and Chanda, 2013). A particular simplicity of such systems is that the maximum number of ligands that can bind the protein molecule (i.e., the protein–ligand stoichiometry) is highly unlikely to change with perturbations (unlike Qmax for voltage-dependent system). Thus, by combining mutant cycle analysis with ligand-binding measurements, the contributions of site-specific interaction energies to ligand dependent activation of proteins can also be assessed.
An analytical comparison of nonadditive energies computed by FMC and GIA in example multistate systems
Consider a voltage-dependent ion channel comprising two subunits. For simplicity we assume that the subunits are identical, such that they transfer an equal amount of gating charge during voltage-dependent activation. For such a system we envision three different gating schemes as shown above. In Scheme A1, the two subunits activate sequentially, following a particular order (say subunit 1 activates before subunit 2). The channel is open only when both subunits are activated. The equilibrium constants of the two steps have the same voltage dependence, q, but differ in their voltage-independent components. In Scheme A2, the two subunits activate independently, but activation of one subunit allosterically modulates the activation of the other, by a factor K2/K1, and the channels are open only when both subunits are activated. As in Scheme A1, both K2 and K1 have the same voltage dependence, q, but differ in their voltage-independent components. Scheme A3 is similar to Scheme A2, but with the central difference that channels are open when any one of the two (or both) subunits are activated (or in other words, the channel is closed only when both subunits are deactivated). For all three schemes, K1 and K2 are written as where i = 1 or 2 and represents the voltage-independent component of the equilibrium constant. Furthermore, for all three schemes, the net free-energy change during voltage-dependent activation (i.e., the chemical component of the free-energy difference between the first and last states of the channel, as indicated by the red arrow) is the same: and so is the total charge transferred during voltage-dependent activation: Qmax = 2q.
What is the PO-V–based Boltzmann metric for the free-energy change in each of the three cases? We use the principle that the Boltzmann slope (obtained through fitting a Boltzmann equation to the PO-V curve) can be approximated to the slope of the Hill-transformed PO-V curve at V1/2 (Yifrach, 2004), i.e.,
Using the above equation, the Boltzmann measure of the free-energy change (i.e., ΔGapp = zappFV1/2) for the three schemes can be derived to be
where ΔGapp-1, ΔGapp-2, ΔGapp-3 are the PO-V–based Boltzmann measures of free-energy change for Schemes A1, A2, and A3, respectively. This derivation illustrates two major points. First, although in all three schemes the total free-energy change associated with full-scale activation of the channel is the same, the Boltzmann measure of free-energy change is different in the three cases. Thus, depending on the nature of the gating scheme, the free-energy change deduced from the PO-V curves will be different. Second, and more importantly, for each of the three schemes, ΔGapp is not related to and by simple logarithmic relations (i.e., unlike ΔGnet). This nonlinear dependence may result in false positives when ΔGapp is used to compute interaction energies as is done in FMC applications.
We thank the members of the Chanda laboratory for their helpful comments and Katherine Baldwin for help with preparing Fig. 7.
This project was supported by funds from the National Institutes of Health (grant R01GM084140), Shaw Scientist Award, and Vilas Research Associate Award to B. Chanda.
The authors declare no competing financial interests.
Kenton J. Swartz served as editor.
S. Chowdhury’s present address is Vollum Institute, Oregon Health and Science University, Portland, OR 97239.