text:
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.

## INTRODUCTION

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.

### Electrophysiology

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).

### Data analysis

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:

$IImax=1exp{zappF(V1/2−V)/RT},$

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

$δΔGapp=F(δzappV1/2)2+(zappδV1/2)2,$

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

$ΔΔGFMC={zappFV1/2}WT+{zappFV1/2}S12−{zappFV1/2}S1−{zappFV1/2}S2,$

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

$δΔΔGFMC=(δΔGapp‖W)2+(δΔGapp‖S12)2+(δΔGapp‖S1)2+(δΔGapp‖S2)2,$

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

$VM=∑i=1n−1(Qi+1−Qi)(Vi+1+Vi)2,$

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

$ΔΔGGIA=QmaxF(VM||W+VM||S12−VM||S1−VM||S2).$

The standard error associated with ΔΔGGIA (δΔΔGGIA) was calculated as

$δΔΔGGIA=QmaxF{δ(VM)WT}2+{δ(VM)S1}2+{δ(VM)S2}2+{δ(VM)S12}2,$

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.

### Simulations

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 $Kimin$ to $Kimax,$ and similarly the voltage-dependent parameters also have a defined range of values $zimin$ to $zimax.$ 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 $Ki,mut=Kimin(Kimax/Kimin)R,$ where R is the random number, whereas the voltage-dependent parameters were assigned as $zi,mut=zimin+R(zimax−zimin)$ (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 $K10:$ [0.005–5,000], $K20:$ [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 $K10:$ 5, $K20:$ 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.

## RESULTS

### 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), $Gint(i),$ can be described by the equation

$Gint(i)=(Gi(X,Y)−Gi(0,Y))−(Gi(X,0)−Gi(0,0))…,$
(1)

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

$Gint(i+1)−Gint(i)={(Gi+1(X,Y)−Gi+1(0,Y))−(Gi+1(X,0)−Gi+1(0,0))}−{(Gi(X,Y)−Gi(0,Y))−(Gi(X,0)−Gi(0,0))},$

which may be rewritten as

$ΔΔGi=(ΔGi→i+1(X,Y)−ΔGi→i+1(0,Y))−(ΔGi→i+1(X,0)−ΔGi→i+1(0,0))=−RTln(Ki(X,Y)Ki(0,0)Ki(0,Y)Ki(X,0))=−RTlnΩi….$
(2)

In Eq. 2, ΔΔGi reflects $Gint(i+1)−Gint(i)$ 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:

$ΔΔGnet=∑i=1nΔΔGi=(ΔG1→n+1(X,Y)−ΔG1→n+1(0,Y))−(ΔG1→n+1(X,0)−ΔG1→n+1(0,0))=−RT∑i=1nlnΩi…,$
(3)

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

$ΔΔGtrue=−4RTln(K1WTK1S12K1S1K1S2)−4RTln(K2WTK2S12K2S1K2S2)−RTln(LWTLS12LS1LS2),$

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

$ΔΔGtrue=−4RTln(JWTJS12JS1JS2)−4RTln(DWTDS12DS1DS2)−RTln(LWTLS12LS1LS2),$

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.

## DISCUSSION

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 ($POmax$) 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 ($−RTlnPOmax$) 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

$ΔΔGGIA=ΔΔ(QmaxFVM)−RTln(POmax,S12.POmax,WPOmax,S1.POmax,S2)…$
(4)

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 $(POmax,S12.POmax,W)/(POmax,S1.POmax,S2)$ becomes larger than ∼25, which would imply that $POmax$ 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), $POmax$ 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 $POmax$ 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 $POmin$ of the channel is significantly large, a correction factor amounting to $RTln(1−POmin)$ 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 $POmin$ > 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.

### Concluding remarks

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.

## APPENDIX

### An analytical comparison of nonadditive energies computed by FMC and GIA in example multistate systems

(Scheme A1)
(Scheme A2)
(Scheme A3)

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 $Ki0exp(qFV/RT),$ where i = 1 or 2 and $Ki0$ 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: $ΔGnet=−RT ln{K10 K20},$ 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.,

$zapp=RTF∂∂Vln(PO1−PO)|V1/2.$

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

$ΔGapp−1=−RT(1+2K202K20+K10+K102+4K10K20).ln(2K10K20K10+K102+4K10K20),$
$ΔGapp−2=−2RT(K20+K10+K102+K10K20K20+2K10+2K102+K10K20).ln(K10K20K10+K102+K10K20),$
$ΔGapp−3=−2RT(K20+K10+K102+K10K20K20+2K10+2K102+K10K20).ln(K10K20K10+K102+K10K20),$

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 $K10$ and $K20$ 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.

## Acknowledgments

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.

## References

References
Ackers
,
G.K.
, and
F.R.
Smith
.
1985
.
Effects of site-specific amino acid modification on protein interactions and biological function
.
Annu. Rev. Biochem.
54
:
597
629
.
Aggarwal
,
S.K.
, and
R.
MacKinnon
.
1996
.
Contribution of the S4 segment to gating charge in the Shaker K+ channel
.
Neuron.
16
:
1169
1177
.
Ahern
,
C.A.
, and
R.
Horn
.
2004
.
Specificity of charge-carrying residues in the voltage sensor of potassium channels
.
J. Gen. Physiol.
123
:
205
216
.
Bezanilla
,
F.
, and
C.A.
Villalba-Galea
.
2013
.
The gating charge should not be estimated by fitting a two-state model to a Q-V curve
.
J. Gen. Physiol.
142
:
575
578
.
Bezanilla
,
F.
,
E.
Perozo
, and
E.
Stefani
.
1994
.
Gating of Shaker K+ channels: II. The components of gating currents and a model of channel activation
.
Biophys. J.
66
:
1011
1021
.
Carter
,
P.J.
,
G.
Winter
,
A.J.
Wilkinson
, and
A.R.
Fersht
.
1984
.
The use of double mutants to detect structural changes in the active site of the tyrosyl-tRNA synthetase (Bacillus stearothermophilus)
.
Cell.
38
:
835
840
.
Chamberlin
,
A.
,
F.
Qiu
,
S.
Rebolledo
,
Y.
Wang
,
S.Y.
Noskov
, and
H.P.
.
2014
.
Hydrophobic plug functions as a gate in voltage-gated proton channels
.
111
:
E273
E282
.
Cheng
,
Y.M.
,
C.M.
Hull
,
C.M.
Niven
,
J.
Qi
,
C.R.
Allard
, and
T.W.
Claydon
.
2013
.
Functional interactions of voltage sensor charges with an S2 hydrophobic plug in hERG channels
.
J. Gen. Physiol.
142
:
289
303
.
Chowdhury
,
S.
, and
B.
Chanda
.
2010
.
Deconstructing thermodynamic parameters of a coupled system from site-specific observables
.
107
:
18856
18861
.
Chowdhury
,
S.
, and
B.
Chanda
.
2012a
.
Estimating the voltage-dependent free energy change of ion channels using the median voltage for activation
.
J. Gen. Physiol.
139
:
3
17
.
Chowdhury
,
S.
, and
B.
Chanda
.
2012b
.
Perspectives on: Conformational coupling in ion channels: Thermodynamics of electromechanical coupling in voltage-gated ion channels
.
J. Gen. Physiol.
140
:
613
623
.
Chowdhury
,
S.
, and
B.
Chanda
.
2013
.
Free-energy relationships in ion channels activated by voltage and ligand
.
J. Gen. Physiol.
141
:
11
28
.
Cordero-Morales
,
J.F.
,
V.
Jogini
,
S.
Chakrapani
, and
E.
Perozo
.
2011
.
A multipoint hydrogen-bond network underlying KcsA C-type inactivation
.
Biophys. J.
100
:
2387
2393
.
DeCaen
,
P.G.
,
V.
Yarov-Yarovoy
,
Y.
Zhao
,
T.
Scheuer
, and
W.A.
Catterall
.
2008
.
Disulfide locking a sodium channel voltage sensor reveals ion pair formation during activation
.
105
:
15142
15147
.
DeCaen
,
P.G.
,
V.
Yarov-Yarovoy
,
E.M.
Sharp
,
T.
Scheuer
, and
W.A.
Catterall
.
2009
.
Sequential formation of ion pairs during activation of a sodium channel voltage sensor
.
106
:
22498
22503
.
DeCaen
,
P.G.
,
V.
Yarov-Yarovoy
,
T.
Scheuer
, and
W.A.
Catterall
.
2011
.
Gating charge interactions with the S1 segment during activation of a Na+ channel voltage sensor
.
108
:
18825
18830
.
Di Cera
,
E.
1998
.
Site-specific thermodynamics: Understanding cooperativity in molecular recognition
.
Chem. Rev.
98
:
1563
1592
.
Ding
,
S.
, and
R.
Horn
.
2002
.
Tail end of the s6 segment: Role in permeation in Shaker potassium channels
.
J. Gen. Physiol.
120
:
87
97
.
Faiman
,
G.A.
, and
A.
Horovitz
.
1996
.
On the choice of reference mutant states in the application of the double-mutant cycle method
.
Protein Eng.
9
:
315
316
.
Gagnon
,
D.G.
, and
F.
Bezanilla
.
2010
.
The contribution of individual subunits to the coupling of the voltage sensor to pore opening in Shaker K channels: effect of ILT mutations in heterotetramers
.
J. Gen. Physiol.
136
:
555
568
.
Gleitsman
,
K.R.
,
J.A.
Shanata
,
S.J.
Frazier
,
H.A.
Lester
, and
D.A.
Dougherty
.
2009
.
Long-range coupling in an allosteric receptor revealed by mutant cycle analysis
.
Biophys. J.
96
:
3168
3178
.
Goodey
,
N.M.
, and
S.J.
Benkovic
.
2008
.
Allosteric regulation and catalysis emerge via a common route
.
Nat. Chem. Biol.
4
:
474
482
.
Gupta
,
S.
, and
A.
Auerbach
.
2011
.
Mapping heat exchange in an allosteric protein
.
Biophys. J.
100
:
904
911
.
Hille
,
B.
2001
.
Ion Channels of Excitable Membranes
. Third edition.
Sinauer Associates
,
Sunderland, MA
.
814 pp
.
Horovitz
,
A.
1996
.
Double-mutant cycles: a powerful tool for analyzing protein structure and function
.
Fold. Des.
1
:
R121
R126
.
Horovitz
,
A.
, and
A.R.
Fersht
.
1990
.
Strategy for analysing the co-operativity of intramolecular interactions in peptides and proteins
.
J. Mol. Biol.
214
:
613
617
.
Horovitz
,
A.
,
E.S.
Bochkareva
,
O.
Yifrach
, and
A.S.
Girshovich
.
1994
.
Prediction of an inter-residue interaction in the chaperonin GroEL from multiple sequence alignment is confirmed by double-mutant cycle analysis
.
J. Mol. Biol.
238
:
133
138
.
Horrigan
,
F.T.
, and
R.W.
Aldrich
.
1999
.
Allosteric voltage gating of potassium channels II. Mslo channel gating charge movement in the absence of Ca2+
.
J. Gen. Physiol.
114
:
305
336
.
Horrigan
,
F.T.
, and
R.W.
Aldrich
.
2002
.
Coupling between voltage sensor activation, Ca2+ binding and channel opening in large conductance (BK) potassium channels
.
J. Gen. Physiol.
120
:
267
305
.
Hoshi
,
T.
,
W.N.
Zagotta
, and
R.W.
Aldrich
.
1994
.
Shaker potassium channel gating. I: Transitions near the open state
.
J. Gen. Physiol.
103
:
249
278
.
Kash
,
T.L.
,
A.
Jenkins
,
J.C.
Kelley
,
J.R.
Trudell
, and
N.L.
Harrison
.
2003
.
Coupling of agonist binding to channel gating in the GABAA receptor
.
Nature.
421
:
272
275
.
Kitaguchi
,
T.
,
M.
Sukhareva
, and
K.J.
Swartz
.
2004
.
Stabilizing the closed S6 gate in the Shaker Kv channel through modification of a hydrophobic seal
.
J. Gen. Physiol.
124
:
319
332
.
Ledwell
,
J.L.
, and
R.W.
Aldrich
.
1999
.
Mutations in the S4 region isolate the final voltage-dependent cooperative step in potassium channel activation
.
J. Gen. Physiol.
113
:
389
414
.
Lee
,
W.Y.
, and
S.M.
Sine
.
2005
.
Principal pathway coupling agonist binding to channel gating in nicotinic receptors
.
Nature.
438
:
243
247
.
Mark
,
A.E.
, and
W.F.
van Gunsteren
.
1994
.
Decomposition of the free energy of a system in terms of specific interactions. Implications for theoretical and experimental studies
.
J. Mol. Biol.
240
:
167
176
.
Miller
,
C.
2012
.
Model-free free energy for voltage-gated channels
.
J. Gen. Physiol.
139
:
1
2
.
Perozo
,
E.
,
R.
MacKinnon
,
F.
Bezanilla
, and
E.
Stefani
.
1993
.
Gating currents from a nonconducting mutant reveal open-closed conformations in Shaker K+ channels
.
Neuron.
11
:
353
358
.
Ranganathan
,
R.
,
J.H.
Lewis
, and
R.
MacKinnon
.
1996
.
Spatial localization of the K+ channel selectivity filter by mutant cycle-based structure analysis
.
Neuron.
16
:
131
139
.
,
E.
, and
O.
Yifrach
.
2007
.
Principles underlying energetic coupling along an allosteric communication trajectory of a voltage-activated K+ channel
.
104
:
19813
19818
.
Schoppa
,
N.E.
, and
F.J.
Sigworth
.
1998
.
Activation of Shaker potassium channels. III. An activation gating model for wild-type and V2 mutant channels
.
J. Gen. Physiol.
111
:
313
342
.
Schoppa
,
N.E.
,
K.
McCormack
,
M.A.
Tanouye
, and
F.J.
Sigworth
.
1992
.
The size of gating charge in wild-type and mutant Shaker potassium channels
.
Science.
255
:
1712
1715
.
Schreiber
,
G.
, and
A.R.
Fersht
.
1995
.
Energetics of protein-protein interactions: analysis of the barnase-barstar interface by single mutations and double mutant cycles
.
J. Mol. Biol.
248
:
478
486
.
Seoh
,
S.A.
,
D.
Sigg
,
D.M.
Papazian
, and
F.
Bezanilla
.
1996
.
Voltage-sensing residues in the S2 and S4 segments of the Shaker K+ channel
.
Neuron.
16
:
1159
1167
.
Serrano
,
L.
,
A.
Horovitz
,
B.
Avron
,
M.
Bycroft
, and
A.R.
Fersht
.
1990
.
Estimating the contribution of engineered surface electrostatic interactions to protein stability by using double-mutant cycles
.
Biochemistry.
29
:
9343
9352
.
Shanata
,
J.A.
,
S.J.
Frazier
,
H.A.
Lester
, and
D.A.
Dougherty
.
2012
.
Using mutant cycle analysis to elucidate long-range functional coupling in allosteric receptors
.
Methods Mol. Biol.
796
:
97
113
.
Sigg
,
D.
2013
.
A linkage analysis toolkit for studying allosteric networks in ion channels
.
J. Gen. Physiol.
141
:
29
60
.
Soler-Llavina
,
G.J.
,
T.H.
Chang
, and
K.J.
Swartz
.
2006
.
Functional interactions at the interface between voltage-sensing and pore domains in the Shaker Kv channel
.
Neuron.
52
:
623
634
.
Wall-Lacelle
,
S.
,
M.I.
Hossain
,
R.
Sauvé
,
R.
Blunck
, and
L.
Parent
.
2011
.
Double mutant cycle analysis identified a critical leucine residue in the IIS4S5 linker for the activation of the CaV2.3 calcium channel
.
J. Biol. Chem.
286
:
27197
27205
.
Wyman
,
J.
, and
S.J.
Gill
.
1990
.
Binding and Linkage: Functional Chemistry of Biological Macromolecules
.
University Science Books
,
Mill Valley, CA
.
330 pp
.
Yang
,
Y.
,
Y.
Yan
, and
F.J.
Sigworth
.
1997
.
How does the W434F mutation block current in Shaker potassium channels?
J. Gen. Physiol.
109
:
779
789
.
Yellen
,
G.
1998
.
The moving parts of voltage-gated ion channels
.
Q. Rev. Biophys.
31
:
239
295
.
Yifrach
,
O.
2004
.
Hill coefficient for estimating the magnitude of cooperativity in gating transitions of voltage-dependent ion channels
.
Biophys. J.
87
:
822
830
.
Yifrach
,
O.
, and
R.
MacKinnon
.
2002
.
Energetics of pore opening in a voltage-gated K+ channel
.
Cell.
111
:
231
239
.
Yifrach
,
O.
,
N.
Zandany
, and
T.
.
2009
.
Examining cooperative gating phenomena in voltage-dependent potassium channels: taking the energetic approach
.
Methods Enzymol.
466
:
179
209
.
Zagotta
,
W.N.
,
T.
Hoshi
, and
R.W.
Aldrich
.
1994a
.
Shaker potassium channel gating. III: Evaluation of kinetic models for activation
.
J. Gen. Physiol.
103
:
321
362
.
Zagotta
,
W.N.
,
T.
Hoshi
,
J.
Dittman
, and
R.W.
Aldrich
.
1994b
.
Shaker potassium channel gating. II: Transitions in the activation pathway
.
J. Gen. Physiol.
103
:
279
319
.
Zandany
,
N.
,
M.
,
I.
Orr
, and
O.
Yifrach
.
2008
.
Direct analysis of cooperativity in multisubunit allosteric proteins
.
105
:
11697
11702
.

Abbreviations used in this paper:

• COVC

cut-open oocyte voltage clamp

•
• FMC

functional mutant cycle

•
• GIA

generalized interaction-energy analysis

•
• TEV

two-electrode voltage

## Author notes

S. Chowdhury’s present address is Vollum Institute, Oregon Health and Science University, Portland, OR 97239.


text:
# Overview

This plugin allows defining macros that process macro argument using external filter program and render its result in Redmine wiki.

For every filter two macros are defined: <macro> and <macro>_include. The first one directly processes its argument using filter, while the second one assumes its argument is wiki page name, so it reads that wiki page content and processes it using the filter.

Macros already bundled with current release are listed below, but adding new ones is typically as easy as adding several lines in plugin config file.

# Installation

1. It’s recommended (but not required) to install popen4 library first as without it plugin is unable to capture stderr output of external command, so it might be hard to debug it if things go wrong.
2. Get sources from github.
3. See Installing a plugin on Redmine site.
4. Copy wiki_external_filter.yml from config folder of plugin directory to the config folder of your redmine installation.
5. After installation it’s strongly recommended to go to plugin settings and configure caching: Administration → Plugins → Wiki External Filter Plugin: Configure and follow instructions. Note that RoR file-based caching suggested by default does not implement proper cache expiration: you should either setup a cron task to clean cache or do it manually from time to time.
6. To successfully use macros with argument expressions, it’s necessary to patch wiki formatting routine so that it preserves macros arguments.

Redmine 1.3.0: apply redmine-1.3.0-macros-escaping.patch to the Redmine core and go to step 7.

Redmine 1.0.2: apply this patch to the Redmine core and go to step 7, there’s no need to patch individual wiki formatters anymore.

Redmine 0.9. if you use Markdown Extra, switch to the following fork of redmine_markdown_extra_formatter that contains all necessary changes. This patch is for older Markdown Extra formatter version and should not be used for new installations.

Note: if you install Markdown Extra formatter - you should enable it in Administration → Settings → General → Text formatting.

If you use wiki formatter other than Markdown Extra (including default Textile formatter) you will have to change your wiki formatter yourselves as follows:

• Change MACROS_RE regexp not to stop too early - in the issue description Textile wiki formatter is mentioned, but in fact this change should be done for whatever wiki formatter you use.
• Change arguments parsing - again, should be done for whatever wiki formatter you use.
• For some of the formatters escaping should be avoided for macro arguments, see how it is done for Markdown Extra wiki formatter.

Preliminary version of patch for Textile support was made available by Yuya Nishihara, according to the author “It works but really messy” but I haven’t tested it so cannot comment on it.

7. To allow passing attachments names as macros arguments Redmine core should be patched accordingly:

Redmine 1.3.0: redmine-1.3.0-attachments.patch.

Redmine 1.0.2: redmine-1.0.2-attachments.patch.

Redmine 0.9. redmine-attachments-in-macros.patch.

Specific filters installation instructions are below.

# Prefedined Macros

## plantuml

PlantUML is a tool to render UML diagrams from their textual representation. It’s assumed that it can be invoked via wrapper /usr/bin/plantuml, here’s its example content:

#!/bin/bash
/usr/bin/java -Djava.io.tmpdir=/var/tmp -jar /usr/share/plantuml/lib/plantuml.jar ${@}  Result is rendered as PNG file. SVG support seems to be under development for PlantUML but so far looks like it’s still unusable. Gentoo ebuild is attached. Example of usage: {{plantuml( Alice -> Bob: Authentication Request alt successful case Bob -> Alice: Authentication Accepted else some kind of failure Bob -> Alice: Authentication Failure opt loop 1000 times Alice -> Bob: DNS Attack end end else Another type of failure Bob -> Alice: Please repeat end )}}  Rendered output: ## graphviz Graphviz is a tool for graph-like structures visualization. It’s assumed that it can be called as /usr/bin/dot. Result is rendered as SVG image or PNG fallback if SVG is not supported by your browser. Example of usage: {{graphviz( digraph finite_state_machine { rankdir=LR; size="8,5" node [shape = doublecircle]; LR_0 LR_3 LR_4 LR_8; node [shape = circle]; LR_0 -> LR_2 [ label = "SS(B)" ]; LR_0 -> LR_1 [ label = "SS(S)" ]; LR_1 -> LR_3 [ label = "S($end)" ];
LR_2 -> LR_6 [ label = "SS(b)" ];
LR_2 -> LR_5 [ label = "SS(a)" ];
LR_2 -> LR_4 [ label = "S(A)" ];
LR_5 -> LR_7 [ label = "S(b)" ];
LR_5 -> LR_5 [ label = "S(a)" ];
LR_6 -> LR_6 [ label = "S(b)" ];
LR_6 -> LR_5 [ label = "S(a)" ];
LR_7 -> LR_8 [ label = "S(b)" ];
LR_7 -> LR_5 [ label = "S(a)" ];
LR_8 -> LR_6 [ label = "S(b)" ];
LR_8 -> LR_5 [ label = "S(a)" ];
}
)}}

Rendered output:

## ritex

Combination of Ritex: a Ruby WebTeX to MathML converter and SVGMath that takes WebTeX formula specification as input and produces SVG file as output.

Both ritex and SVGMath require some patches/wrappers.

Additionally working installation of xmllint from libxml2 with configured MathML catalog is required: for Gentoo use this ebuild.

Gentoo ebuilds for ritex and svgmath are attached.

Example of usage:

{{ritex(
G(y) = \left\{\array{ 1 - e^{-\lambda x} & \text{ if } y \geq 0 \\ 0 & \text{ if } y < 0 }\right.
)}}

Rendered output:

## video and video_url

These macros use ffmpeg to convert any supported video file to FLV format and display it on wiki using FlowPlayer flash player. video macro takes file path on server as its input, as well as attachments names from current wiki page, while video_url expects full URL to the video to convert & show.

Splash images for videos are generated automatically from the first frame of the video.

Multiple videos per page are supported, player instance is attached to the selected video as in this example.

Required packages installed:

• ffmpeg
• RMagick
• wget - for video_url only.

Example of usage:

{{video(video_attachment_name.avi)}}

or

{{video(/full/path/to/the/file/on/the/server/video.avi)}}

or

{{video_url(http://www.example.com/video.avi)}}

Rendered output (before player is embedded by clicking on the image, using Flowplayer video demo file):

## fortune

Fortune is a simple program that displays a random message from a database of quotations.

Not strictly a filter on its own (as it does not require any input), but it plays nice with external filtering approach and is fun to use, hence it’s here

Example of usage:

{{fortune}}

Rendered output:

# Third-Party macros

## rail

This macro provides Syntax diagrams support, thanks to Siamak Haschemi, see his blog for details and macro code.

# Writing new macros

New macros can easily be added via wiki_external_filter.yml config file.

Every macro may have multiple commands processing the same input - for example for video macro two commands are used: first one extracts thumbnail and second one converts the video.

Commands use standard Unix approach for filtering: input is fed to the command via stdin and output is read on stdout. If command return status is zero, content type is assumed to be of content_type specified in config, otherwise it’s assumed to be plain error text together with stderr content.

You can use prolog/epilog config parameters to add standard text before/after actual macro content passed to filter.

Additionally, cache_seconds parameter specifies the number of seconds commands output result should be cached, use zero to disable caching for this macro.

The way filter output is visualized is controlled via app/views/wiki_external_filter/macro_*.html.erb files. The view to use is selected by template macro option in config. The view can use all commands outputs for particular macro.

replace_attachments tells plugin that it should parse the text passed to the macro and replace all occurrences of strings matching attachments names with their physical paths on disk.

Macro argument is de-escaped via CGI.unescapeHTML call prior to being fed to filter.

# Current bugs/issues

1. Either Redmine core (if you use default wiki engine) or your custom wiki engine plugin requires patching to get things work. In fact, the whole wiki formatting design as of now seems to be quite messy.
2. SVG support is more complex it should have been if all browsers had played by the rules - currently quite some trickery with different XHTML elements/CSS tricks is used to show SVGs properly in major browsers. Of course, there’s not that much that can be done for IE as it does not support SVG at all, but now at least the plugin substitutes raster fall-back image for IE if it is available.
3. For formula support, theoretically ritex alone is sufficient if you have MathML-capable browser, however in practice there are too many issues with this approach: for example Firefox (actually the onlt MathML-capable browser so far, it seems) requires specific DOCTYPE additions that Redmine currently lacks; additionally, Redmine emits text/html, while Firefox expects text/xml in order to parse MathML. Changing content type alone is not sufficient as Redmine HTML output does not pass more strict checks required for XML output. Hence, the double conversion (WebTeX to MathML and then MathML to SVG) is necessary. Once (if ever?) MathML support matures in other browser, possibly this can be revisited.
4. SVGs could have been embedded into HTML page directly (thus allowing to use redmine links there) but I’m afraid there are similar problems as with MathML embedding attempts.
5. RoR caching support is a mess: no way to expire old files from file-based cache??? Are you joking???

format: plain


text:
# 1 Installation

To install the package, please use the following.

if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("mbkmeans")

# 2 Introduction

This vignette provides an introductory example on how to work with the mbkmeans package, which contains an implementation of the mini-batch k-means algorithm proposed in (Sculley 2010) for large single-cell sequencing data. This algorithm runs the _k_means iterations on small, random subsamples of data (mini batches’’) rather than the entire dataset. This both speeds up computation and reduces the amount of data that must be contained in memory.

The main function to be used by the users is mbkmeans. This is implemented as an S4 generic and methods are implemented for matrix, Matrix, HDF5Matrix, DelayedMatrix, SummarizedExperiment, and SingleCellExperiment.

Most of this work was inspired by the MiniBatchKmeans() function implemented in the ClusterR R package and we re-use many of the C++ functions implemented there.

Our main contribution here is to provide an interface to the DelayedArray and HDF5Array packages so that mbkmeans algorithm only brings into memory at any time the “mini-batches” needed for the current iteration, and never requires loading the entire dataset into memory. This allows the user to run the mini-batch k-means algorithm on very large datasets that do not fit entirely in memory. mbkmeans also runs seamlessly on in-memory matrices provided by the user, and doing so significantly increases the speed of the algorithm compared to the standard k-means algorithm ((Hicks et al. 2020)). A complete comparison of mbkmeans with other implementations can be found in (Hicks et al. 2020).

The motivation for this work is the clustering of large single-cell RNA-sequencing (scRNA-seq) datasets, and hence the main focus is on Bioconductor’s SingleCellExperimentand SummarizedExperiment data container. For this reason, mbkmeans assumes a data representation typical of genomic data, in which genes (variables) are in the rows and cells (observations) are in the column. This is contrary to most other statistical applications, and notably to the stats::kmeans() and ClusterR::MiniBatchKmeans() functions that assume observations in rows.

We provide a lower-level mini_batch() function that expects observations in rows and is expected to be a direct replacement of ClusterR::MiniBatchKmeans() for on-disk data representations such as HDF5 files. The rest of this document shows the typical use case through the mbkmeans() interface; users interested in the mini_batch() function should refer to its manual page.

## 2.1 Example dataset

To illustrate a typical use case, we use the pbmc4k dataset of the TENxPBMCData package. This dataset contains a set of about 4,000 cells from peripheral blood from a healthy donor and is expected to contain many types or clusters of cell.

We would note that mbkmeans is designed for very large datasets, in particular for datasets large enough that either the full data matrix can’t be held in memory, or the computations necessary cannot be computed. This vignette works with a small dataset to allow the user to quickly follow along and understand the commands. (Hicks et al. 2020) shows the example of using mbkmeans with a dataset with 1.3 million cells, which is a more appropriate size for observing improved memory usage and improved speed.

First, we load the needed packages.

library(TENxPBMCData)
library(scater)
library(SingleCellExperiment)
library(mbkmeans)
library(DelayedMatrixStats)

Note that in this vignette, we do not aim at identifying biologically meaningful clusters (that would entail a more sophisticated normalization of data and dimensionality reduction), but instead we only aim to show how to run mini-batch k-means on a HDF5-backed matrix.

We normalize the data simply by scaling for the total number of counts using scater and select the 1,000 most variable genes and a random set of 100 cells to speed-up computations.

tenx_pbmc4k <- TENxPBMCData(dataset = "pbmc4k")

set.seed(1034)
idx <- sample(seq_len(NCOL(tenx_pbmc4k)), 100)
sce <- tenx_pbmc4k[, idx]

#normalization
sce <- logNormCounts(sce)

vars <- rowVars(logcounts(sce))
names(vars) <- rownames(sce)
vars <- sort(vars, decreasing = TRUE)

sce1000 <- sce[names(vars)[1:1000],]

sce1000
## class: SingleCellExperiment
## dim: 1000 100
## assays(2): counts logcounts
## rownames(1000): ENSG00000090382 ENSG00000204287 ... ENSG00000117748
##   ENSG00000135968
## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
## colnames: NULL
## colData names(12): Sample Barcode ... Date_published sizeFactor
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

# 3mbkmeans

The main function, mbkmeans(), returns a list object including centroids, WCSS_per_cluster (where WCSS stands for within-cluster-sum-of-squares), best_initialization, iters_per_initiazation and Clusters.

It takes any matrix-like object as input, such as SummarizedExperiment, SingleCellExperiment, matrix, DelayedMatrix and HDF5Matrix.

In this example, the input is a SingleCellExperiment object.

res <- mbkmeans(sce1000, clusters = 5,
reduceMethod = NA,
whichAssay = "logcounts")

The number of clusters (such as k in the k-means algorithm) is set through the clusters argument. In this case, we set clusters = 5 for no particular reason. For SingleCellExperiment objects, the function provides the reduceMethod and whichAssay arguments. The reduceMethod argument should specify the dimensionality reduction slot to use for the clustering, and the default is “PCA”. Note that this does not perform PCA but only looks at a slot called “PCA” already stored in the object. Alternatively, one can specify whichAssay as the assay to use as input to mini-batch k-means. This is used only when reduceMethod option is NA. See ?mbkmeans for more details.

## 3.1 Choice of arguments

There are additional arguements in mbkmeans() that make the function more flexible and suitable for more situations.

### 3.1.1 Batch size

The size of the mini batches is set through the batch_size argument. This argument gives the size of the samples of data brought into memory at any one time, and thus constrains the memory usage of the alogrithm.

The blocksize() function can be used to set the batch size to the maximum allowed by the available memory. It considers both the number of columns in the dataset and the amount of RAM on the current matchine to calculate as large of a batch size as reasonable for the RAM available to the session. The calculation uses get_ram() function in benchmarkme package. See the benchmarkme vignette for more details.

batchsize <- blocksize(sce1000)
batchsize
## [1] 100

In this simple example, because the whole data fits in memory, the default batch size would be a single batch of size 100.

(Hicks et al. 2020) provides a comparison of the effect of batch size for large datasets (i.e. where the full data cannot be stored in memory). Their results show the algorithm to be robust to a large range of batch sizes (1,000-10,000 cells): accuracy of the results were equivalent if the batch size was above 500-1000 cells, with no noticable advantage in memory usage or speed as long as the batches were smaller than roughly 10,000 cells (corresponding to a subsample of the full data matrix whose memory foot print is not noticeable on most modern computers). For this reason, we would err on the larger range of batch size suggested in (Hicks et al. 2020) – 10,000 cells. Furthermore, since the initialization uses by default the same number of cells as batch size, this gives a robust sample for determining the initial start point.

We would note that it is better to make choices on batch size based on choosing the absolute number of cells in a batch rather than choosing batch size based on a percentage of the overall sample size. If batch size is chosen based on the percentage of cells, very large datasets will use very large batches of cells, negating the memory and speed gains. Moreover such large batches are not needed for good performance of the algorithm. The continual resampling of the batches works like stochastic gradient descent and allows for (local) minimization of the objective function (Sculley 2010) with relatively small batch sizes, as demonstrated on single-cell data in the work of (Hicks et al. 2020). If anything, the size of the batch should be governed more by the complexity of the problem (e.g. the number of underlying true clusters or subtypes), rather than the number of overall cells in the dataset.

### 3.1.2 Initialization

The performance of mini-batch k-means greatly depends on the process of initialization. We implemented two different initialization methods:

1. Random initialization, as in regular k-means;
2. kmeans++, as proposed in (Arthur and Vassilvitskii 2007). The default is “kmeans++”.

The percentage of data to use for the initialization centroids is set through the init_fraction argument, which should be larger than 0 and less than 1. The default value is given so that the proportion matches the value of batch_size, converted to a proportion of the data.

res_random <- mbkmeans(sce1000, clusters = 5,
reduceMethod = NA,
whichAssay = "logcounts",
initializer = "random")
table(res$Clusters, res_random$Clusters)
##
##      1  2  3  4  5
##   1  0  0  0  0  1
##   2  0  0  1  0  0
##   3  4  0 23  0  0
##   4  0  0  0  0  4
##   5  0  5  0 14 48

Note that large values of init_fraction will result in that proportion of the data being held in memory and used in the initialization of the kmeans++ algorithm. Therefore, large values of init_fraction will require a great deal of memory – potentially more than that used by the actual clustering algorithm. Indeed keeping the fraction used constant for larger and larger datasets will result in commiserate numbers of cells being used in the initialization and increasing the memory usage (similar to the issues discussed above for batch size). Therefore, we recommend that users keep the default and make changes to the batch_size parameter as desired, ensuring consistent memory usage across both the initialization stage and the actual clustering stage of the algorithm. As mentioned in the discussion above on batch_size, for this reason we would recommend batch_size on the larger range suggested as reasonable in (Hicks et al. 2020) so as to improve the initialization.

## 3.2 Running mbkmeans with multiple values of $$k$$

The main parameter to set in $$k$$-means and its variants is the number of clusters $$k$$. mbkmeans is quick enough to rerun the clustering algorithm with different number of clusters even in very large datasets.

Here, we apply mbkmeans with $$k$$ from 5 to 15 and select the number of clusters using the elbow method (i.e., the value that corresponds to the point of inflection on the curve). We note that this is just a rule of thumb for selecting the number of clusters, and many different methods exist to decide the appropriate value of $$k$$.

To speed up the computations, we will cluster on the top 20 PCs.

sce1000 <- runPCA(sce1000, ncomponents=20)

ks <- seq(5, 15)
res <- lapply(ks, function(k) {
mbkmeans(sce1000, clusters = k,
reduceMethod = "PCA",
calc_wcss = TRUE, num_init=10)
})

wcss <- sapply(res, function(x) sum(x$WCSS_per_cluster)) plot(ks, wcss, type = "b") From the plot, it seems that k = 12 could be a reasonable value. # 4 Comparison with k-means Note that if we set init_fraction = 1, initializer = "random", and batch_size = ncol(x), we recover the classic k-means algorithm. res_full <- mbkmeans(sce1000, clusters = 5, reduceMethod = NA, whichAssay = "logcounts", initializer = "random", batch_size = ncol(sce1000)) res_classic <- kmeans(t(logcounts(sce1000)), centers = 5) table(res_full$Clusters, res_classic$cluster) ## ## 1 2 3 4 5 ## 1 0 7 0 7 0 ## 2 0 0 14 43 0 ## 3 2 0 0 0 0 ## 4 16 0 0 0 10 ## 5 0 0 0 1 0 Note however, that since the two algorithms start from different random initializations, they will not always converge to the same solution. Furthremore, if the algorithm uses batches smaller than the full dataset (batch_size<ncol(x)), then mbkmeans will not converge to the same result as kmeans even with the same initialization. This is due to the fact that kmeans and mbkmeans have different search paths for an optimum and the problem is non-convex. For a comparison of memory usage, speed, and accuracy of mbkmeans with kmeans for the non-trivial case where the algorithm uses batches, see (Hicks et al. 2020), where it is demonstrated that mbkmeans provides significant speed improvements and uses less memory than kmeans, without loss of accuracy. # 5 Use with bluster The bluster package provides a flexible and extensible framework for clustering in Bioconductor packages/workflows. The clusterRows() function controls dispatch to different clustering algorithms. In this case, we will use the mini-batch k-means algorithm to cluster cells into cell populations based on their principal component analysis (PCA) coordinates. mat <- reducedDim(sce1000, "PCA") dim(mat) ## [1] 100 20 In the following three scenarios, we pass the number of cluster centers as a number using the centers=5 argument. In addition, we can also pass other mbkmeans arguments (e.g. batch_size=10). In the first two cases, we only return the cluster labels. However, in the third scenario, we ask for the full mbkmeans output. library(bluster) clusterRows(mat, MbkmeansParam(centers=5)) ## [1] 2 5 5 3 5 4 3 1 1 3 5 2 5 1 1 1 3 5 1 3 2 1 5 5 1 3 1 1 5 1 3 5 4 1 2 1 4 ## [38] 5 5 4 5 4 5 4 5 3 1 3 5 5 5 2 1 3 2 5 4 1 4 3 5 4 1 1 3 5 2 4 5 4 1 5 4 1 ## [75] 5 3 1 5 3 1 2 1 4 1 2 5 4 5 1 1 5 1 4 2 4 5 2 5 1 2 ## Levels: 1 2 3 4 5 clusterRows(mat, MbkmeansParam(centers=5, batch_size=10)) ## [1] 3 3 3 3 3 3 3 5 5 3 3 1 3 5 5 5 3 3 5 3 4 5 3 3 5 3 5 5 3 5 3 3 4 5 1 5 4 ## [38] 3 3 3 3 3 3 3 3 3 5 3 3 3 3 1 5 3 1 3 3 5 4 3 3 4 5 5 3 3 3 3 3 3 2 3 3 5 ## [75] 3 3 5 3 3 5 4 5 3 5 1 3 4 3 5 5 3 5 3 3 3 3 3 3 5 4 ## Levels: 1 2 3 4 5 clusterRows(mat, MbkmeansParam(centers=5, batch_size=10), full = TRUE) ##$clusters
##   [1] 5 5 5 5 5 5 5 1 1 5 5 2 5 1 1 1 5 5 3 5 5 1 5 5 1 5 1 1 5 1 5 5 5 1 2 1 5
##  [38] 5 5 5 5 5 5 5 5 5 1 5 5 5 5 2 1 5 2 5 5 1 5 5 5 5 1 1 5 5 5 5 5 5 1 5 5 1
##  [75] 5 5 1 5 5 1 5 1 5 1 2 5 5 5 1 1 5 1 5 5 5 5 4 5 1 5
## Levels: 1 2 3 4 5
##
## $objects ##$objects$mbkmeans ##$objects$mbkmeans$centroids
##            [,1]       [,2]      [,3]       [,4]      [,5]       [,6]       [,7]
## [1,] -15.010417 -0.7962079  4.968180 -2.7679840 -0.605855 -0.2194566 -1.4958129
## [2,]   1.398788 15.2268393 -4.914318 -8.7031508 -0.624902 -2.6562225 -0.8620272
## [3,] -17.543300 -0.3406707  6.138011 -4.3145411  1.079512  0.6031462  3.6588221
## [4,]   3.505305  9.0542211 -4.008002 -0.8464335 -1.494348 -5.7764435 -3.4572454
## [5,]   4.953126 -5.2662492 -1.885661  0.1000841  1.303848 -0.6517419 -0.3327930
##            [,8]       [,9]      [,10]      [,11]     [,12]      [,13]
## [1,] -0.6686381  0.7806568 -1.1359341  1.1528491 -1.407013  1.5457278
## [2,] -3.2219319 -3.0456175 -0.2216553 -1.8928053  5.403923 -2.9907297
## [3,]  8.2619290 -4.0679099 -3.2909212  5.7171778  8.625866  4.8241775
## [4,] -3.3455593 -2.0061982 -1.3108509 -9.4799142  1.971986  9.3333830
## [5,] -0.4513868 -1.7769361 -1.2359026  0.4449523  0.151892 -0.9781383
##           [,14]      [,15]      [,16]       [,17]     [,18]     [,19]
## [1,]  1.6981873  3.2334236 -0.7828917  0.11174810 0.7921813  2.095090
## [2,] -1.0326219 -0.3752647  2.2624371  2.48013343 2.6617969 -2.687306
## [3,]  0.2712364 -1.0289170 -5.9773702 -2.34843263 1.5557200 -3.698316
## [4,]  0.8668148  1.5455810 -1.8839526 -4.15320464 2.9828210  2.413375
## [5,]  2.0870510 -1.4041611 -0.3987755 -0.06157393 0.6216171 -0.435778
##            [,20]
## [1,] -0.12527754
## [2,] -2.21443762
## [3,]  2.67774700
## [4,]  0.06607878
## [5,] -0.77234827
##
## $objects$mbkmeans$WCSS_per_cluster ## numeric(0) ## ##$objects$mbkmeans$best_initialization
## [1] 1
##
## $objects$mbkmeans$iters_per_initialization ## [,1] ## [1,] 17 ## ##$objects$mbkmeans$Clusters
##   [1] 5 5 5 5 5 5 5 1 1 5 5 2 5 1 1 1 5 5 3 5 5 1 5 5 1 5 1 1 5 1 5 5 5 1 2 1 5
##  [38] 5 5 5 5 5 5 5 5 5 1 5 5 5 5 2 1 5 2 5 5 1 5 5 5 5 1 1 5 5 5 5 5 5 1 5 5 1
##  [75] 5 5 1 5 5 1 5 1 5 1 2 5 5 5 1 1 5 1 5 5 5 5 4 5 1 5

While many of the times clustering will be performed on in-memory data (for example after PCA), the clusterRows function also accepts any matrix-like object (in-memory or on-disk data representation) that mbkmeans accepts.

For example, we can cluster using mbkmeans on the logcounts assay in the sce1000 object (similar to above), which is a DelayedMatrix object from the DelayedArray package.

Note: We transpose the matrix as the clusterRows function expects observations along the rows and variables along the columns.

logcounts(sce1000)
## <1000 x 100> sparse matrix of class DelayedMatrix and type "double":
##                      [,1]      [,2]      [,3] ...     [,99]    [,100]
## ENSG00000090382  0.000000  1.482992  0.000000   . 5.4330604 0.0000000
## ENSG00000204287  2.198815  0.000000  0.000000   . 5.3387971 0.0000000
## ENSG00000163220  1.762993  1.482992  0.000000   . 2.1176311 2.0642508
## ENSG00000143546  0.000000  1.135438  1.113099   . 0.3825221 0.0000000
## ENSG00000019582  3.402730  1.762861  1.733874   . 6.2009459 2.0642508
##             ...         .         .         .   .         .         .
## ENSG00000145779 0.0000000 0.6766523 0.0000000   . 0.9342195 0.0000000
## ENSG00000106948 0.0000000 1.1354382 0.0000000   . 0.9342195 1.3735556
## ENSG00000155660 0.0000000 0.6766523 0.0000000   . 0.3825221 1.3735556
## ENSG00000117748 0.0000000 0.0000000 0.0000000   . 0.0000000 0.0000000
## ENSG00000135968 0.0000000 0.0000000 0.0000000   . 0.3825221 0.0000000
clusterRows(t(logcounts(sce1000)), MbkmeansParam(centers=4))
##   [1] 2 2 2 1 2 2 1 3 3 1 2 2 2 3 3 3 1 2 3 1 2 3 2 2 3 1 3 3 2 3 1 2 2 3 2 3 2
##  [38] 2 2 2 2 2 2 2 2 1 3 1 2 2 2 2 3 1 2 4 2 3 2 1 2 2 3 3 1 2 2 2 2 2 3 2 2 3
##  [75] 2 1 3 2 1 3 2 3 2 3 2 2 2 2 3 3 2 3 2 2 2 2 2 2 3 2
## Levels: 1 2 3 4

# 6 Working with on-disk data in Bioconductor

This vignette is focused on the use of the mbkmeans package. For the users interested in learning more about working with HDF5 files and on-disk data in Bioconductor, we recommend the excellent DelayedArray workshop, presented by Peter Hickey at Bioc2020.

# References

Arthur, David, and Sergei Vassilvitskii. 2007. “K-Means++: The Advantages of Careful Seeding.” In Proceedings of the Eighteenth Annual Acm-Siam Symposium on Discrete Algorithms, 1027–35. Society for Industrial; Applied Mathematics.

Hicks, Stephanie C., Ruoxi Liu, Yuwei Ni, Elizabeth Purdom, and Davide Risso. 2020. “Mbkmeans: Fast Clustering for Single Cell Data Using Mini-Batch K-Means.” bioRxiv. https://doi.org/10.1101/2020.05.27.119438.

Sculley, David. 2010. “Web-Scale K-Means Clustering.” In Proceedings of the 19th International Conference on World Wide Web, 1177–8. ACM.


text:
If xy=z and (x+33)(y-29)=z, then what x and y? : GMAT Quantitative Section
Check GMAT Club Decision Tracker for the Latest School Decision Releases https://gmatclub.com/AppTrack

 It is currently 26 Feb 2017, 06:11

### GMAT Club Daily Prep

#### Thank you for using the timer - this advanced tool can estimate your performance and suggest more practice questions. We have subscribed you to Daily Prep Questions via email.

Customized
for You

we will pick new questions that match your level based on your Timer History

Track

every week, we’ll send you an estimated GMAT score based on your performance

Practice
Pays

we will pick new questions that match your level based on your Timer History

# Events & Promotions

###### Events & Promotions in June
Open Detailed Calendar

# If xy=z and (x+33)(y-29)=z, then what x and y?

Author Message
TAGS:

### Hide Tags

Intern
Joined: 29 Mar 2013
Posts: 3
Followers: 0

Kudos [?]: 0 [0], given: 0

If xy=z and (x+33)(y-29)=z, then what x and y? [#permalink]

### Show Tags

29 Mar 2013, 22:46
The solution of the above escapes me, as I know the answer, but I would like to know howg to solve it.

Regards Jan
Intern
Joined: 29 Mar 2013
Posts: 3
Followers: 0

Kudos [?]: 0 [0], given: 0

Re: If xy=z and (x+33)(y-29)=z, then what x and y? [#permalink]

### Show Tags

30 Mar 2013, 01:24
The 'z' is of no conseque, but there should be two solutions, one being x and y are two real numbers.
Intern
Joined: 29 Mar 2013
Posts: 3
Followers: 0

Kudos [?]: 0 [0], given: 0

Re: If xy=z and (x+33)(y-29)=z, then what x and y? [#permalink]

### Show Tags

30 Mar 2013, 04:21
Sorry for my ambiguity. Z is just the same product of these equations. I have been told, that there to solutions, one of which are two real numbers e.g 5 and 8.
Senior Manager
Status: Prevent and prepare. Not repent and repair!!
Joined: 13 Feb 2010
Posts: 274
Location: India
Concentration: Technology, General Management
GPA: 3.75
WE: Sales (Telecommunications)
Followers: 9

Kudos [?]: 90 [0], given: 282

Re: If xy=z and (x+33)(y-29)=z, then what x and y? [#permalink]

### Show Tags

01 Apr 2013, 11:52
this is like bull in a china shop.. running everywhere? not gmat like!
_________________

I've failed over and over and over again in my life and that is why I succeed--Michael Jordan
Kudos drives a person to better himself every single time. So Pls give it generously
Wont give up till i hit a 700+

Intern
Joined: 15 Jan 2013
Posts: 39
Concentration: Finance, Operations
GPA: 4
Followers: 0

Kudos [?]: 28 [0], given: 6

Re: If xy=z and (x+33)(y-29)=z, then what x and y? [#permalink]

### Show Tags

02 Apr 2013, 05:18
jjwho wrote:
The solution of the above escapes me, as I know the answer, but I would like to know howg to solve it.

Regards Jan

$$xy = (x + 33)(y - 29)$$ = > $$33y - 29x = 29.33$$. It is a line and will have infinite solutions.
Re: If xy=z and (x+33)(y-29)=z, then what x and y?   [#permalink] 02 Apr 2013, 05:18
Similar topics Replies Last post
Similar
Topics:
DS question: What is the value of x^2 - y^2? 0 14 Jul 2016, 13:24
3 x,y,z are even positive integers. Does it mean all are different? 3 13 Dec 2015, 05:32
1 If x, y, and z are consecutive integers, is x+y+z divisble by 3? 4 22 Feb 2015, 03:49
Clarity on what " x,y,z are consecutive integers" means 3 16 Sep 2011, 10:19
23 If integers x, y and z are greater than 1 what is the value 23 04 Aug 2008, 10:05
Display posts from previous: Sort by


text:
topos theory

# Contents

## Idea

An abelian sheaf is a sheaf with values in an abelian category which usually is itself, or is taken to be embedded in, a category of complexes in an abelian category.

In light of the Dold-Kan correspondence this means that abelian sheaves can usefully be regarded as special cases of simplicial presheaves and in particular the corresponding derived category of abelian sheaves, traditionally mainly investigated in terms of sheaf cohomology, is analogous to the homotopy category of abelian infinity-stacks. In this way, via Dold-Kan, plain abelian sheaves already go a long way towards (abelian) Higher Topos Theory, which is one way of understanding the relevance of the concept of abelian sheaves.

For instance Deligne cohomology, which classifies higher abelian gerbes (certain infinity-stacks) with connection), is the sheaf cohomology of a certain class of sheaves with values in abelian complexes. This is understood conceptually by realizing that after embedding complexes of abelian sheaves – via Dold-Kan – into general simplicial sheaves, a complex of abelian sheaves becomes an abelian $\infty$-prestack and the computation of its sheaf cohomology corresponds to passing to its infinity-stackification.

## Properties

### Projective objects

See at projective object the section Existence of enough projectives.

## References

A basic textbook introduction begins for instance around Definition 1.5.6 of

A detailed textbook discussion is in section 18 of

category: sheaf theory

Revised on November 24, 2013 07:58:16 by Urs Schreiber (89.204.137.79)