A thermodynamic approach to studying allosterically regulated ion channels such as the large-conductance voltage- and Ca2+-dependent (BK) channel is presented, drawing from principles originally introduced to describe linkage phenomena in hemoglobin. In this paper, linkage between a principal channel component and secondary elements is derived from a four-state thermodynamic cycle. One set of parallel legs in the cycle describes the “work function,” or the free energy required to activate the principal component. The second are “lever operations” activating linked elements. The experimental embodiment of this linkage cycle is a plot of work function versus secondary force, whose asymptotes are a function of the parameters (displacements and interaction energies) of an allosteric network. Two essential work functions play a role in evaluating data from voltage-clamp experiments. The first is the conductance Hill energy WH[g], which is a “local” work function for pore activation, and is defined as kT times the Hill transform of the conductance (G-V) curve. The second is the electrical capacitance energy WC[q], representing “global” gating charge displacement, and is equal to the product of total gating charge per channel times the first moment (VM) of normalized capacitance (slope of Q-V curve). Plots of WH[g] and WC[q] versus voltage and Ca2+ potential can be used to measure thermodynamic parameters in a model-independent fashion of the core gating constituents (pore, voltage-sensor, and Ca2+-binding domain) of BK channel. The method is easily generalized for use in studying other allosterically regulated ion channels. The feasibility of performing linkage analysis from patch-clamp data were explored by simulating gating and ionic currents of a 17-particle model BK channel in response to a slow voltage ramp, which yielded interaction energies deviating from their given values in the range of 1.3 to 7.2%.

Ion channels are pore-forming proteins that regulate (gate) the flux of selected ions (K+, Na+, Ca2+, Cl) across the cell membrane in response to a variety of external forces (environmental stimuli) that include membrane potential (V), ligand chemical potential (μ), temperature (T), membrane tension, and light energy (Hille, 2001; Nagel et al., 2002). In many channels, gating is strongly influenced by a ring of accessory domains acting as efficient conveyers of thermal, mechanical, or chemical energy to the pore’s central gating apparatus. An essential source of mechanical energy for many ion channels is Q-V work. The product of Q (membrane-specific charge) and V (membrane potential) is one of several canonical force-displacement terms in the system energy. Using voltage-clamp techniques, it is possible to study charge movement in the channel regulatory apparatus. This can be done by measuring Q directly, typically in the form of a “gating” current Ig = dQ/dt, which, when integrated, contributes a data point in the equilibrium Q-V curve. An indirect method of studying voltage gating and other means of environmental sensitivity is to measure the conductance G, a specific marker of pore activation. Equilibrium curves other than the standard G-V and Q-V curves can be useful, for example G-T and G-μ (for a list of commonly used parameters that appear in this paper, please see Table 1).

The main objective when analyzing equilibrium data are to identify elements that respond to external forces, and to quantify their interactions with the catalytic unit, which in ion channels is the conducting pore. Quaternary descriptions of enzyme function have been useful in studying regulatory proteins—most notably hemoglobin, widely considered the poster child of protein allosteric theory. Modeling hemoglobin using sophisticated variants of the classical Monod–Wyman–Changeux (MWC) equation (Monod et al., 1965) has achieved impressive insight into its allosteric machinery (Eaton et al., 2007). A K+ channel whose regulation, at a basic level, is formulaically similar to that of hemoglobin (though mechanistically distinct), is the large-conductance voltage- and Ca2+-dependent (BK) channel. The BK channel derives its voltage dependence from four voltage-sensing (J) domains located within the membrane electric field, and also to a small degree from the pore (L) itself. Calcium sensors (K) are situated beneath the membrane plane, where they form a structure known as the gating ring. With voltage and Ca2+ sensors each outnumbering the pore 4:1, a significant mechanical leverage can be exerted on the pore gate. Combining this with the observation that pore opening is effectively (though perhaps not strictly) a binary process, one is lead to the concerted (MWC) model as a natural choice for allosteric regulation of BK.

The relatively large number of regulatory domains (four voltage sensors and possibly eight or more Ca2+-binding sites) generates a multitude of possible channel states, the bulk of which appear to be readily accessible through experimental application of V and [Ca2+]i. A landmark study by Horrigan and Aldrich (2002) firmly established the allosteric nature of the BK channel through a wide range of kinetic and equilibrium measurements, culminating in a model with nine interacting gating domains or particles (one pore, four voltage sensors, and four Ca2+-binding sites). Since then, derivatives of the Horrigan–Aldrich (HA) model containing an even greater number of gating particles have been proposed to account for various experimental findings (Zeng et al., 2005; Qian et al., 2006; Horrigan and Ma, 2008; Sweet and Cox, 2008; Pantazis et al., 2010; Savalli et al., 2012).

This paper describes thermodynamic methods that can be used to characterize allosteric networks in ion channels. The theory proposed here has firm roots in the pioneering work by Jeffries Wyman (1964) on ligand-binding allosteric proteins. Wyman referred to interactions between liganded “protomers” or “particles” (regulatory units of proteins) as “linkage,” and argued on thermodynamic grounds that linkage must be reciprocal; that is, the energy WXY linking the activated state of a protomer Y to that of a neighboring protomer X equals the energy of the reverse interaction WYX. Linkage theory has generally focused on interactions between ligand-binding protomers, which is reflected in a notational style that overlooks the existence of other environmental forces. Here, Wyman’s original theory is generalized for use with any external force through the means of a mathematical device referred to as the “particle potential” η, which simply is the free energy of particle activation expressed as a sum of force-displacement canonical pairs, related to the particle’s equilibrium constant by K = exp(−η/kT). An important use of η is to derive the equilibrium curve of particle activation, which equals ∂Φ/∂η, where Φ is the system energy. Armed with the partition function Z (a statistically weighted sum of channel states expressed as a polynomial function of equilibrium constants and allosteric factors—the equivalent of Wyman’s binding polynomial), linkage relations between regulatory protomers are easily derived. This is a powerful and unifying concept that provides a framework for understanding the thermodynamics of multimodal, allosterically regulated proteins such as BK.

A key variable in the current treatment is the “work function,” an energy specifically assigned to the principal component (A) of a linkage relation. The principal component may be a canonical displacement (for example, total gating charge Q), in which case we speak of a “global” variable, or it could be a marker of protomer activation (for example, the pore open state), in which case we are dealing with a “local” variable (Di Cera, 1990). The work function WA is formally defined as the change in system energy Φ required to reversibly activate A. It is as if one possessed an imaginary lever for this purpose and could measure the effort WA = ΦA(+) − ΦA(−) required to very slowly change A from the resting to the activated state. One can also envision other levers whose function is to activate secondary components (B) that might interact with A. These secondary levers do not require a measurement of effort because the goal is to determine the interaction energy WAB and not the total energy of activating both components (= WA + WB + WAB). This is helpful because, in many instances, it is not possible to measure WB, but it is possible to activate components B through a lever operation. The lever operation can in principle be performed by any channel-modifying agent, for example, a mutation or a toxin. Here, however, we will consider only generalized forces FB such as voltage or calcium chemical potential as secondary agents, as these can be smoothly and quantitatively controlled in the laboratory.

Using a combination of lever pulls, the energy of interaction WAB between components A and B can be obtained by twice measuring WA (i.e., pulling the A lever), once with component B in the resting state (WA,B(−)) and a second time with component B in the activated state (WA,B(+)). The difference between these two energies yields WAB, the logic being that if A and B are allosterically linked, activating B will assist (or retard) activation A by the amount WAB. This form of linkage analysis can be concisely summarized by the expression WAB = BΔWA = BΔAΔΦ, where the “lever operator” Δ() is defined (using component A as an a example) by AΔ() = ()A(+) − ()A(−). In plain English, a lever operation takes the difference of the operand (the quantity in parentheses) evaluated at the extreme limits of a second quantity. The process of lever pulls is neatly represented by a four-state thermodynamic cycle, in which each corner represents a limiting energy state in the equation WAB = (ΦA(+)B(+) + ΦA(−)B(−) − ΦA(+)B(−) − ΦA(−)B(+)). In both pathways starting in the resting state ΦA(−)B(−) and ending in the fully activated state ΦA(+)B(+), one leg constitutes the work function WA, whereas the other leg is defined by the lever operation BΔ (see Fig. 2 B). The reciprocity principle demands that the alternative cycle, which uses WB and AΔ as its legs, should yield the same interaction energy; in other words, WAB = WBA. Here we confine the discussion to pairwise interactions. More complicated cycles with eight or more states (see Sadovsky and Yifrach, 2007), such as those that describe the interaction WABC (in which the activation state of a third component C influences WAB), are not covered here but represent a straightforward extension of the theory.

The allosteric relationship between components A and B is established by plotting WA against FB, referred to here as a linkage plot. The typical appearance of a linkage plot is shown in Fig. 1. The work function WA rises steeply in the region of B activation, with a displacement equal to WAB. If the principal component A is itself not sensitive to FB, the plot plateaus at the two extremes. Otherwise, the extreme behavior is characterized by sloping asymptotic lines whose derivative m is the generalized displacement of A linked to FB (for example, if FB = V, then m = ΔqA, the gating charge of A). There is a strong resemblance between the linkage plot and the Hill plot used to study cooperative binding in hemoglobin and other allosteric proteins. This is not accidental as we shall see, although the methods presented here focus on the limiting asymptotes of the plot rather than on the steeper middle portion, the slope of the latter used primarily as a measure of cooperativity related to the Hill coefficient (Hill, 1910; Yifrach, 2004).

To move the discussion beyond the pulling of imaginary levers to the more practical realm, we must specify work functions that can be measured in the laboratory. In the BK channel, and in other voltage-dependent ion channels studied under voltage clamp, two such functions stand out as being essential, based on their connections to the equilibrium curves G-V and Q-V. The first of these relates to pore opening, for which G is a local marker. The work function for pore activation is, by definition, WL = ΦL(+) − ΦL(−). It will be shown later that WL is derived from the G-V curve through WL = −kTln[G/(GmaxG)]. Given the similarity of this formula to the traditional Hill plot of ligand binding (Wyman, 1964; Eaton et al., 2007), and also its relationship to G, −WL will be referred to as the “conductance Hill energy” WH[g]. Defining WH[g] as the negative value of −WL generates positive asymptotic slopes characteristic of the traditional Hill plot. Other markers such as fluorescent labels can be used to define additional local variable Hill functions, for example, WJ and WK, which are the work functions for voltage and Ca2+ sensor activation, respectively. The function lnε used in “χ-analysis” of mutational effects on cooperativity (Chowdhury and Chanda, 2010) is precisely equal to WH[g] apart from a factor of kT.

The second essential work function is a global variable obtained from the Q-V curve. It is the reversible work Wq of moving a channel’s gating charge q = Q/N, where Q is total measured charge, and N is the number of channels. Wq is derived by integrating Vdq over the saturating range of q given by [0, ΔqT]. Changing the integrand from q to V yields Wq = ΔqTVM, where VM = ∫VfqdV is the mean value of V with respect to the capacitance distribution fq = (dq/dV)/ΔqT. Because the work function for q depends in part on the global “capacity” of the channel to store charge, it is referred to here as the “electrical capacitance” energy WC[q], which we define as WC[q] = −Wq, again to conform to the positive sloping convention of the traditional Hill plot. The capacitance energy of other thermodynamic displacements can be similarly defined. For example, Wn = ΔnTμM describes the reversible energy of ligand binding, where μM is the mean value of the chemical potential μ integrated over the “binding capacity” dn/ (Di Cera et al., 1988), and n is the ligand number.

The summed energy Wqn of interactions between voltage-sensing and Ca2+-sensing elements in the BK channel is equal to the expression μΔWq = ΔqT μΔVM, whose evaluation in practice involves measuring the shift in VM under limiting conditions of “zero” and “saturating” calcium concentration, denoted by μ(−) and μ(+). Later, it will be demonstrated how a slow voltage ramp can be used to elicit the capacitance Cq = dQ/dV, from which VM can be obtained. Wyman’s reciprocity principle assures us that VΔWn = μΔWq, or more succinctly, that Wqn = Wnq. In other words, measuring the mean ligand activity μM as a function of V could in principle be used to confirm the value of Wqn obtained from the μ dependence of VM, provided one is able to measure Ca2+-binding curves.

Linkage plots of global work functions can generate sloping asymptotes resembling those of Hill plots, provided that at least one of the contributors to the capacitance is activated by both principal and secondary forces. An example would be a Ca2+-binding site located a small electrical distance d inside the membrane potential. If the resulting charge movement 2d can be detected as a component of the Q-V curve, plotting WC[q] versus μ would demonstrate asymptotic slopes equal to Δn, the number of Ca2+ bound per site. However, in most models of BK, there is the assumption of strict modality separation between voltage and Ca2+ sensors, so that the extremes of the linkage plot are horizontal plateaus. The absence of sloping asymptotes is a useful feature if one wants to economize on data collection, as only two measurements (WC[q],μ[−] and WC[q],μ[+]) are required to measure Wqn = μΔWq.

With the two essential work functions, WH[g] and WC[q], one is able to obtain the three core interaction energies linking the L, J, and K protomers of the BK channel. This will be explicitly demonstrated for the well-known nine-particle HA model (Scheme 2 in this paper), but is also true for more complex models of BK, where J and K protomers may contain a multitude of internal interactions (Scheme 3). The results from Monte Carlo simulation using patch-clamp conditions for an explicit model containing 17 particles (Scheme 4) will be analyzed to help determine the feasibility of performing linkage analysis on BK channels in the laboratory. A slow voltage-ramp protocol was used for the simulated experiments, with the advantage over conventional step-pulse protocols that only single sweeps are needed to generate “quasi-equilibrium” Q-V or G-V curves.

Several issues arise when considering linkage analysis in ion channels, a fair number of which will be discussed in this paper. For example, how do local work functions behave when there are multiple copies of the particle of interest, particularly if the copies self-interact? This will be addressed when we consider a multiple interacting five-particle model (Scheme 1), and again for a more complex 20-state particle (Scheme 5) in which binary conductance is theorized to be the result of strong pore subunit interactions. A somewhat speculative scenario is when coupling energies themselves are functions of external forces. This leads to distortions in the linkage plots that will be explored in the context of the 17-particle model. A third question speaks to the appropriateness of thermodynamic methods for channels other than BK. Is linkage analysis appropriate for use in all K+ channels? The answer is probably not. Members of the Kv class or voltage-dependent K+ channels are likely not amenable to linkage analysis of the voltage sensor–pore interaction, presumably because large interaction energies reduce the accessibility of most open-state configurations, a statement that is compatible with the current view that the open state of the Shaker K+ channel can be visited only if all four voltage sensors are activated (Schoppa and Sigworth, 1998; Ledwell and Aldrich, 1999; Horn et al., 2000; Gagnon and Bezanilla, 2009). However, it is likely that other allosterically regulated channels possess a more accessible configurational space. An interesting example that will be briefly explored is the temperature sensitivity of a class of transient receptor potential (TRP) channels.

### Data simulations

Q-V and G-V curves for the 17-particle model (Scheme 3) were computed from particle activation curves (Eqs. 39ac) using the modeling program Berkeley Madonna. Monte Carlo patch-clamp simulations of Scheme 3 kinetics were performed using custom software written in C. Each Monte Carlo experiment consisted of 100 traces of summed currents from n = 1,000 ± 50 channels as well as contributions from thermal noise and membrane-based currents (leak and capacitance). At any time t in the simulation, the state of the channel was represented by a bitwise array containing the activation states (0 for resting and 1 for activated) of the 17 regulatory particles. All gating particles were in their resting state at the beginning of the voltage protocol (t = 0). After allowing the system to equilibrate during a 50-ms constant-voltage prepulse at −400 mV, the voltage was slowly ramped with a speed m = 1 mV/ms to an endpoint of 300 mV, ending with a 50-ms constant-voltage postpulse also at 300 mV. Intervals τ between transition events were randomly determined by assigning each particle i a waiting time di based on its state-dependent transition rate αi. During constant-voltage segments of the simulation, or if the particle charge Δqi = 0, di was obtained from the formula −lnrn/αi, where rn is a uniform random number drawn between 0 and 1 obtained with the long-period (>2 × 1018) random number generator of L’Ecuyer with Bays-Durham shuffle (Press et al., 1992). During the ramp phase, the di of charged particles was obtained from the formula (2kT/mΔqi)ln[1−(mΔqilnrn)/(2kTαoi)] (Sigg and Bezanilla, 1997), which reflects the nonstationary value of the forward and backward particle rate constants αi and βi, which are of the form αoi exp(ΔqiV/2kT) and βoi exp(−ΔqiV/2kT), respectively. The method for determining αi and βi are given in the next section. The particle with the shortest di was allowed to transition, after which the simulation time was advanced by di. If between t and t + di there was a change in pulse protocol occurring at t1, the simulation time reverted to t1, and the process was resumed for the new set of conditions. The gating current for each transition was computed as the filtered impulse response (cutoff frequency fc) using a Gaussian filter and discretized (sampling frequency fs = 5fc) in such a way that the discrete integral of the response equaled the gating charge displacement Δqi (Sigg et. al., 1999). Because transition intervals were stored as floating-point numbers, and discretization was performed in concert with filtering, multiple transitions could occur within the same sampling time increment 1/fs without introducing error. The unitary ionic current i was generated as a sequence of filtered unitary impulses (alternating between areas −1 and +1) connected to each closing and opening event, which was then integrated into a telegraph signal and multiplied by the pore conductance gL and driving force (VVeq). The value of the equilibrium potential Veq in ramp simulations was 50.0 ± 2.5 mV. Gaussian-distributed white random noise was added to each trace of ∼1,000 channels. The noise amplitude was adjusted so that, after digital filtering, the mean-square value of the resulting current (Inoise) was equal to the Nyquist value 4kTBg, where B is the low-pass bandwidth of a Gaussian filter (B ≈ 1.064 fc; see Crouzy and Sigworth, 1993), and the conductance g was made up of contributions from a constant “leak” conductance gleak, assumed to be 50.0 ± 2.5 pS (≈ 20 Gohm−1), and a fluctuating term NOgL, where NO was the number of open channels at a given time. Each simulated current trace (gating and ionic) contained the following terms:

$I=Nig+Ni+Inoise+mCpatch+gleakV,$

where Cpatch is the patch capacitance (= 3.00 ± 0.15 fC), and the ionic current i was set to zero when simulating gating currents. Perfect series resistance compensation was assumed. The parameters N, Veq, gleak, and Cpatch were randomly generated to within a range of ±5% of their stated values using a uniform random deviate. This was done at the beginning of each experiment. The randomized values of these patch-specific parameters were not known at the time of analysis. In addition to ramp simulations, conventional square-pulse protocols were used to simulate charge-per-channel experiments in which gating and ionic currents were generated for a 10-ms test (ON) pulse to 300 mV from a resting potential Vr = −300 mV, ending with a tail (OFF) pulse to 0 mV.

### Kinetics of the 17-particle model

Although the emphasis of this paper is on equilibrium methods, kinetic parameters for the 17-particle model were required to perform ramp simulations. The chosen scheme was consistent with known macro-kinetics of the channel while preserving equilibrium properties. Time constants of decay τ for both ionic and gating currents in BK were obtained from Horrigan and Aldrich (2002). In their study, maximum value τmax of the voltage-sensor decay time, estimated from the τ versus V plot of the rapid component of gating current, was ∼5 ms, whereas ionic current kinetics, reflecting pore activation, demonstrated a maximum τmax in the decay time of ∼0.05 ms. Both voltage-sensor particles in Scheme 3 were assigned the same decay time. Ca2+ binding was assumed to have the same τmax as voltage-sensor activation. The condition of maximum τ = (α + β)−1 was obtained when forward (α) and backward (β) particle activation rates were both equal to (2τmax)−1. To satisfy equilibrium and kinetic requirements, particle rate constant were computed according to αi = (2τmax)−1(Zi(+)/Zi(−))1/2 and βi = (2τmax)−1(Zi(−)/Zi(+))1/2, where Zi(−) and Zi(+) are the limiting channel partition functions corresponding to the resting and activated states of the ith particle, obtained during ramp simulation from a look-up table continuously updated to account for changes in V.

### Data analysis

Simulated experiments consisted of 100 traces each of ionic and gating currents from n = 1,000 ± 50 channels, taken at [Ca2+]i = 10−3 µM and 103 µM. Estimates for gleak and Cpatch were obtained by analyzing short segments (2–40 ms) of current before (I(−)) and after (I(+)) the beginning of the test pulse (t1). For a square pulse of height ΔV, gleak = (I(+)I(−))/ΔV, where the current was averaged over each segment. In ramp experiments, initial estimates of gleak and Cpatch were simultaneously obtained by fitting I to the following function: I (t < t1) = I(−); I (tt1) = I(−) + mCpatch + gleak(V(t) − Vr). Fits were obtained using the solver application in Excel (Microsoft). The values of gleak and Cpatch, in addition to Veq (in ionic current experiments), were further refined through manual adjustment and fitting of asymptotes to achieve satisfactory shapes for both the Q-V curve = ∫(Nig/m)dV and the conductance Hill plot WH[g] = kTln(G/(GmaxG), where G = I/(VVeq). To reduce the deleterious effects of Nyquist current noise on Q-V and WH[g] plots, post-filtering of currents was performed after initial correction for gleak and Cpatch with a symmetric Gaussian filter to an adjusted cutoff frequency of 0.1 kHz. After adjusting baselines for the ramp-generated gating currents as above, the value of VM was estimated from the following integral (see Eq. 29):

$VM=∫VigdV∫igdV.$

Integration was performed using Riemann sums, which was adequate for the sampling rate used (δV = 0.2 mV). In WH[g] versus V plots, straight lines were fit to the positive and negative asymptotes, and energies of interaction were obtained from the height differences. This was done for both [Ca2+]i = 10−3 and 103 µM. The pore gating charge ΔqL was obtained from the slope of the lower (negative) V asymptote in the high Ca2+ plot. In “charge-per-channel” experiments, the time-dependent variance Var(I) of the ionic current was plotted against I and fitted to the equation Var(I) = iII2/N (Sigworth, 1980) for both test and tail currents, with N averaged over the two estimates. The total gating charge ΔqT was obtained from the integrated positive gating current ΔQ = ∫NigdV divided by the measured n value.

Parameter estimates in the form of mean ± SEM (SEM) were obtained from 10 experiments in each set of conditions. Two-tailed Student’s t tests were used to assign p-values to the deviation of estimates from their expected values. Differences with p-values of <0.05 were considered statistically significant.

### Online supplemental material

The partition function and linkage relationships for Scheme 3 (see Eqs. 37ad) are derived. The use of linkage analysis to discrimination between mechanisms of phenotype reversal in chimeric temperature-sensing TRP channels is illustrated in Fig. S1. Supplemental text and Figs. S1 and S2 are available at http://www.jgp.org/cgi/content/full/jgp.201210859/DC1.

### Thermodynamics

We consider a patch-clamp experiment measuring N channels at fixed temperature and pressure. The bath and membrane contain several buffered or reservoir species (divalent ions, H+, water, lipids) and additional species found in fixed numbers (the channel, other ancillary proteins, unbuffered ligands). The fundamental equation of thermodynamics for this system is:

$dU=TdS–PdV+VdQ+ΣbμbdNb+ΣfμfdNf,$
(1)

where U is energy, T is temperature, S is entropy, P is pressure, V (extensive) is volume, V (intensive) is voltage, Q is charge, and μb and μf are the chemical potentials of buffered and fixed-numbered species present in numbers Nb and Nf. There may be other mechanical terms not included in this equation, such as those related to membrane tension, which find use in the study of stretch-activated channels (Markin and Sachs, 2004). The chemical species are additionally divided by what side of the membrane (intracellular or extracellular) they inhabit.

We are interested in a modified Gibbs energy, obtained by applying Legendre transforms to all controlled variables, including buffered species, namely, G[V,{μb}] = UST + PVQV − ΣbμbNb (notation adapted from Alberty, 2002). Integrating Eq. 1 at constant values of the intensive variables, we obtain G[V,{μb}] = Σfμf Nf, indicating that the free energy is a function of fixed-numbered species only. Assuming that the only fixed-numbered species of interest is the BK channel, we obtain G[V,{μb}] = BK. Following Wyman (1975), we equate the system energy Φ with μBK. Writing the Gibbs–Duhem equation for Eq. 1 as −SdT + VdPQdV − ΣbNbbNdΦ = 0 and solving for dΦ, we obtain:

$dΦ= –sdT+vdP–qdV–Σbnbdμb,$
(2)

where previously extensive variables are now written in lower case to indicate that they have been divided by N and made intensive. As a consequence, Eq. 2 contains only intensive variables. Assuming constant T and P, and allowing only one of the buffered chemical species to vary (internal Ca2+ in the case of BK), we simplify Eq. 9 to derive the system equation for BK, placing angle brackets around displacement variables q and n to indicate averaging:

$dΦ=−〈q〉dV−〈n〉dμ.$
(3)

In Eq. 3, the Ca2+ chemical potential μ = μο + kTlnaCa is a function of internal calcium activity aCa (in practice, assumed to be [Ca2+]i, the molar concentration) in units of micromoles/liter (µM). The standard chemical potential μο can be made zero by choosing the reference concentration to be 1.0 µM.

We assume that linkage relations between canonical pairs q-V and n-μ are entirely mediated by the channel protein, or if external sources of linkage exist, such as surface charge effects (Moczydlowski et al., 1985), they can be corrected. This assumption is formalized by use of the “bridge” equation,

$Φ=−kTlnZ,$
(4)

where Z is the channel-specific partition function. Here, kT has its usual thermodynamic significance (= 25.3 meV at 20°C).1 The partition function, or weighted sum of states, is expressed in terms of the channel’s gating particles as follows:

$Z=∑j,k,l=0nj,nk,nlΘj,k,lexp(−(jηJ+kηK+lηL)kT),$
(5)

where the weights Θ are functions of combinatorial and allosteric factors, and integers j, k, and l are the activation numbers for protomers J, K, and L, ranging from 0 to a maximum of 4 (nJ, nK) or 1 (nL). The variables pertaining to the voltage sensor (= J) and gating ring (= K) are written in bold type to suggest a conformational complexity requiring composite variables, whereas the pore (L) for now has only one activated state. To keep the equations simple, there will be an implied summation when expressions of energies and displacements of composite (bolded) variables are considered (for example, ∂Φ/∂ηJ = ∂Φ/∂ηJ1 + ∂Φ/∂ηJ2 +…), and a product with composite allosteric factors (for example J = J1J2…).

A basic premise of allosteric theory is that protomers behave as quasi-independent units, with a given protomer X possessing a “system” energy ΦX that is defined as the energy landscape of activation experienced by X when all other protomers are at rest. Defining the “particle potential” ηX as the change in ΦX required to activate the X protomer (ηXXΔΦX), and evaluating ηX at constant T, P, V (voltage), and μ, we obtain:

$ηX=ΔGX−ΔqXV−ΔnXμ,$
(6)

where ΔqX and ΔnX are charge and Ca2+ displacements of activation, and ΔGX = Δ(U/N)X − ΔsXT + ΔvXP is the Gibbs energy of particle activation (not to be confused with the macroscopic conductance G and an allosteric factor by the same name). The specific form of Eq. 6 in models of BK can be derived from a “truth table” (Table 2), which lists properties of individual particles. For example, ηL = ΔGX − ΔqLV, because the pore is assumed to be purely voltage dependent.

An “accessible” allosteric system is defined as one for which the weights Θj,k,l in Eq. 5 are nonzero for any combination of j, k, and l. The partition function Z can be converted into polynomial form:

$Z=∑j,k,l=0nj,nk,nlΘj,k,lJjKkLl,$
(7)

where the particle equilibrium constants J, K, and L are related to their respective particle potentials through “local bridge” equations (compare with Eq. 4) ηJ = −kTlnJ, ηK = −kTlnK, and ηL = −kTlnL. From these relations and the distribution of states provided by Eq. 7, a statistical mechanical version of the system equation (Eq. 3) can be derived:

$kTdlnZ=−〈j〉dηK−〈j〉dηJ−〈l〉dηL,$
(8)

where 〈j〉 = ∂Φ/∂ηJ = ∂lnZ/∂lnJ is the mean number of activated J particles, and the other activation numbers are similarly given by 〈k〉 = ∂Φ/∂ηK = ∂lnZ/∂lnK and 〈l〉 = ∂Φ/∂ηL = ∂lnZ/∂lnL. The activation numbers establish a link between local and global descriptions of the system. For example, we can derive the mean gating charge displacement 〈q〉 = −(∂Φ/∂V)μ as a function of particle transition charges by applying the chain rule:

$〈q〉=−∂Φ∂ηJ(∂ηJ∂V)μ−∂Φ∂ηK(∂ηK∂V)μ−∂Φ∂ηL(∂ηL∂V)μ.$
(9)

Evaluating the derivatives on the right side with the help of Eq. 6 yields:

$〈q〉=〈j〉ΔqJ+〈k〉ΔqK+〈l〉ΔqL.$
(10)

An analogous expression, the mean number of bound calcium ions, is:

$〈n〉=〈j〉ΔnJ+〈k〉ΔnK+〈l〉ΔnL.$
(11)

The activation curves 〈j〉, 〈k〉, and 〈l〉 are easily derived from the polynomial form of the partition function, Eq. 7.

### Defining the terms “protomer” and “particle”

The words “protomer” and “particle,” referring to the basic building blocks of allosteric networks, have been used almost interchangeably so far in this paper. However, it is useful to propose a distinction between the two terms to clarify the role of complex domains with multiple activation states. The word “particle” will hereafter refer to a single activated protomer state. In the case of a two-state protomer, such as the pore with partition function Z = 1 + L, the terms “protomer” and “particle” are synonymous. However, a composite protomeric unit possesses more than one activated state. This may be because it undergoes successive transitions during activation, as is thought to be true for the S4 segment in Kv channels (Schoppa and Sigworth, 1998; Gamal El-Din et al., 2010; Henrion, et al., 2012), in which case the partition function is written Z = 1 + J1 + J2 +…; or it is composed of multiple interacting two-state particles, in which case the partition function is (for two particles) Z = (1 + J1) + J2(1 + J1G), with G as a coupling factor. In both cases, the particles are connected to the equilibrium constants J1, J2, etc., and can exist either in the real sense (representing a physical structure or binding site) or virtual sense (a state in a series of activated states). Note that for very large G and small J2, the second class of partition function has the same form as the first. Thus, the sequentially activating protomer can be mathematically viewed as a special case of the interacting particle protomer. Regardless of the form of the partition function, particle activation numbers can always be derived from Z using 〈j1〉 = ∂lnZ/∂lnJ1, 〈j2〉 = ∂lnZ/∂lnJ2, etc.

Just as one (composite) protomer can be made up of multiple particles, the converse may also be true, in which multiple protomers interact so strongly as to activate in concert, thus representing a single particle. Such a scenario explains binary conductance in a K+ channel pore constructed from four identical P subunits. The corresponding partition function is Z = 1 + P4 from which the classical Hill equation (Hill, 1910) can be derived as H = (1/4)∂lnZ/∂lnP = P4/(1 + P4). The two-state partition function Z = 1 + L generates the Boltzmann equation, B = ∂lnZ/∂lnL = L/(1 + L), from which we conclude that in the concerted activating pore, the equilibrium constant is L = P4.

Incorporating a composite protomer into an allosteric scheme is implemented by enabling its constituent particles, whether real or virtual, to interact with other particles. Each configurational state becomes a term within the partition function of Eq. 7, from which equilibrium curves, such as those described by Eqs. 10 and 11, can be derived. In this paper, we will only consider networks composed of simple two-state protomers because this is traditionally what is used to model BK.

### Contracted partition functions (CPFs)

The concept of limiting states, in which the channel adopts a constrained configuration or sub-scheme, either through the use of a mathematical device such as Hill transformation or with the application of an external force, was introduced earlier in the discussion of thermodynamic cycles, and has been instrumental in the understanding of thermodynamic models of BK in relation to experimental data (Rothberg and Magleby, 1999; Horrigan and Aldrich, 2002). To illustrate the method, we examine pore activation, which when visualized in single-channel recordings, appears as a binary stochastic sequence resembling a telegraph signal, with abrupt transitions occurring between closed (C) and open (O) states. The partition function for this two-state process is given by:

$Z = ZC4+LZO4.$
(12)

The first and second terms on the right side of Eq. 12 are “limiting” partition functions (LPFs) for the closed and open sub-schemes of the channel, respectively, and can be symbolized by ZL(−) and ZL(+). Generally speaking, an LPF can be reduced by dividing through common factors to a “contracted” partition function (CPF; see Di Cera, 1990). CPFs differ from LPFs in that their energetic point of reference is not necessarily the least activated state of the channel. The CPFs in Eq. 12 are ZC and ZO, which represent subunit configuration states consistent with the closed and open states of the pore. They are raised to the fourth power in Eq. 12 because we assume for now that subunits do not self-interact.

A useful feature of CPFs is that suitably chosen ratios yield expectation values (distribution-based averages) of coupling factors. To demonstrate this, consider the following truncated model of BK consisting of a single (voltage-dependent) pore particle (L) interacting with four Ca2+-sensing K particles (Scheme 1). The equilibrium scheme for this five-particle model—which, incidentally, is mathematically equivalent to the original MWC model (Monod et al., 1965; Cox et al., 1997)—is, using “linkage notation” (see Table 3 and Appendix):

(SCHEME 1)

In Scheme 1, pore opening (L) simultaneously increases all four equilibrium constants K by a factor of C. The partition function for Scheme 1 is given by Eq. 12, where ZC = (1 + K), and ZO = (1 + KC). The ratio of the two CPFs yields the expectation value of C derived from the distribution of K particle states:

$ZOZC=(1+KC1+K)=〈C〉K.$
(13)

The Ca2+ dependence in Eq. 13 arises from the equilibrium constant K = exp(−ηK/kT) = exp[(−ΔGK + ΔnKμ)/kT]. Recalling that [Ca2+]i = exp(μ/kT), this can be rearranged to read K = ([Ca2+]i/Kd)ΔnK, where Kd = exp(ΔGKnKkT) is the intrinsic disassociation constant of the K-binding site. If a single ion is bound during activation, then ΔnK = 1, and we obtain the usual expression K = [Ca2+]i/Kd. Because K grows monotonically with increasing μ, it is correct to substitute 〈Cμ for 〈CK on the right side of Eq. 13. The value of 〈Cμ changes smoothly from 1 to C when varying μ over the saturating range [μ(−), μ(+)]. Switching to energy units, we conclude that kTln〈Cμ varies from 0 to −ΔWC.

The preceding example is representative of the correspondence between the sigmoidal component of a linkage plot and the logarithm of a CPF ratio. The ability to measure a linkage plot in practice is therefore contingent upon finding a work function that contains the desired CPF ratio. This is the subject of the remainder of the Results section, where work functions based on the G-V and Q-V curves are used to analyze linkage in a series of increasingly complex models of BK.

### The conductance Hill energy WH[g]

The Hill energy is defined as the work function of a single particle. It was stated in the Introduction that the pore-specific Hill energy WL can be derived from the G-V curve. This will now be proved. Starting with the definition WL = ΦL(+) − ΦL(−), we invoke the bridge equation (Eq. 4) to obtain WL = −kTlnZL(+)/ZL(−). As discussed in the preceding section, the partition function can be written as the sum of two LPFs: Z = ZL(−) + ZL(+). Recognizing that each term in Z is a state probability after normalizing with Z, we write the probability PO of channel opening as G/Gmax = ZL(+)/Z. Combining the above equations, and recalling that the conductance Hill energy was defined as WH[g] = −WL, we obtain, as advertised, the following Hill equation:

$WH[g]= kTln(GGmax−G).$
(14)

In channels possessing a binary pore and noninteracting subunits, Z is given by Eq. 12, from which we derive the following statistical mechanical version of Eq. 14 (see also Chowdhury and Chanda, 2010):

$WH[g]= 4kTln(ZOZC)−ηL.$
(15)

From the discussion in the previous section, we know that WH[g] + ηL varies in sigmoid fashion from 0 to −4ΔWC in the range [μ(−), μ(+)]. Applying the lever operator μΔ to Eq. 15, we therefore obtain the following practical formula for the total interaction energy between L and four K particles, seen in the energy Hill plot as the rise in the sigmoidal component:

$Δμ(WH[g]+ηL)=−4ΔWC.$
(16)

It should be noted that, because the pore in Scheme 1 is not intrinsically Ca2+ dependent, its particle potential ηL = ΔGL − ΔqLV is independent of μ; i.e., μΔηL = 0. Under these circumstances, ηL can be left out of Eq. 16.

Plotting WH[g] versus V instead of μ yields a straight line given by WH[g] = −ηL. This is the Hill plot of the pore particle in the absence of interactions with other voltage sensors. Under conditions of “zero” Ca2+, symbolized by μ(−), the line has slope m = ΔqL and V-intercept VL,μ(−) = ΔGLqL, consistent with ηL = −ΔqL(VVL,μ(−)). If we increase Ca2+ to saturating levels, then m is unchanged but the V-intercept VL,μ(+) = (ΔGL + 4ΔWC)/ΔqL is left-shifted in the setting of positive allosteric interaction (ΔWC < 0). This single-particle behavior is the source of sloping asymptotes in Hill plots.

### The ligand-binding Hill plot for Scheme 1

The conductance Hill plot (Eq. 14) was derived from the equilibrium curve (G-V) of a single particle L. Suppose we could measure Ca2+ binding in Scheme 1 through an assay in which the signal B increments by b each time a Ca2+ ion is bound, and we would like to construct an energy Hill plot around the four K particles. Would it be correct to use the following analogue of Eq. 14:

$WH[b]=kTln(BBmax−B),$
(17)

where B = Nbk〉 (maximum value: Bmax = 4Nb)? The answer is yes, supported by the everyday use of conventional Hill plots in ligand-binding assays. To demonstrate this explicitly in the case of Scheme 1, one could evaluate 〈k〉 = ∂lnZ/∂lnK using Eq. 12, but it is more instructive to parse Z differently by expanding it in powers of K, using the following alternative version of the Scheme 1 linkage diagram:

$(L)→K(LC)→K2(LC2)→K3(LC3)→K4(LC4)$

(SCHEME 1, alternative form)

The corresponding partition function (equivalent to Eq. 12) is given by:

$Z = Z(0) + 4Z(1)K + 6Z(2)K2 + 4Z(3)K3 + Z(4)K4,$
(18)

where the weights Z(k) = (1 + LC k) are CPFs corresponding to integer steps in K activation (k = 0…4). Defining ZK = ∂Z/∂K, we derive from Eq. 17 the following:

$WH[b]=kTln(KZK4Z−ZK).$
(19)

Evaluating Eq. 19 with respect to Eq. 18, we obtain the four-particle analogue to Eq. 15:

$WH[b] = kTln(Z(1)+3Z(2)K+3Z(3)K2+Z(4)K3Z(0)+3Z(1)K+3Z(2)K2+Z(3)K3)−ηK.$
(20)

It is possible to regain from Eq. 20 the simple binary ratio of CPFs that made it possible to extract the value of −4ΔWC from Eq. 15, by noting that the log quantity in parentheses converges to binary CPF pairs for either K→0 or K→∞ (corresponding to the limiting conditions μ(−) and μ(+)), yielding:

$WH[b]μ(−)= kTln(Z(1)Z(0))−ηK,$
(21a)
$WH[b]μ(+)= kTln(Z(4)Z(3))−ηK.$
(21b)

We can best understand Eqs. 21a and 21b by recognizing that under limiting conditions, only one K particle (the first or last in the sequence) is statistically able to activate at a time, thus recapitulating the single-particle Hill plot. Applying the voltage lever operation VΔ to either equation yields the anticipated value of −ΔWC (however, notice the absence of the factor 4 compared with Eq. 16). Because the particle potential ηK = ΔGK − ΔnKμ is a linear function of μ, the plot of WH[b] versus μ demonstrates sloping asymptotes (m = ΔnK), as illustrated in Fig. 1 (substituting WH[b] for WA, and μ for FB). Assuming ΔnK = 1, the expressions for the μ intercepts of the two asymptotic lines are μK,μ(−) = ΔGK and μK,μ(+) = ΔGK + ΔWC, from which the disassociation constant Kd is obtained through Kd = exp(μK/kT). We see that, stemming from linkage with the pore, there are two values for Kd (Xia et al., 2002), one for the un-liganded (free) state and the other for the bound state (Kd value left-shifted for the bound state if ΔWC < 0).

### Interacting K particles in Scheme 1

Thus far, we have assumed that each K particle is unaware of the state of the other three K particles except as inferred from the activation state of the pore. In reality, there is extensive intersubunit contact between Ca2+-binding domains that make up the BK gating ring (Yuan et al., 2010). Despite evidence (Qian et al., 2006) that intrasubunit interactions are more relevant than intersubunit interactions for the regulation of pore gating, it is useful to consider the existence of self-interacting K particles, a scenario that can be described by the “square” variant of the Koshland, Némethy, and Filmer mechanism of allosterism (Koshland et al., 1966). We postulate a modification of Eq. 18, in which the CPFs Z(k) are multiplied by successive products of interaction factors Bk. As a result, the original Z(k) are changed to B0Z(0), B0B1Z(1), B0B1B2Z(2), and so forth. The zeroth factor B0 can be set to unity without loss of generality. Substituting these modified CPFs into Eq. 21, we see that the Koshland, Némethy, and Filmer scheme adds energies of interaction kTlnB1 and kTlnB4 to Eqs. 21a and 21b, respectively. Therefore, in evaluating μΔWH[b] for Scheme 1, we obtain—in addition to the K–L interaction energy −ΔWC and the divergent term −μΔηK = ΔnK(μ(+)μ(−))—a third, self-interaction term kTln(B4/B1), representing the difference in the work needed to activate the first and last K particles.

### Summary of energy Hill plot analysis

Summarizing the results so far, we conclude that energy Hill analysis yields three distinct energies: a diverging term (η) related to particle displacement, a heterologous term comprised of the interaction energies between the principal and secondary particles, and a homologous term describing self-interaction among copies of the principal particle. The first of these terms can always be distinguished from the other two through its asymptotic behavior, whereas the heterologous and homologous interaction terms combine to generate the sigmoidal component of the Hill plot (Fig. 1). Note that, in global linkage analysis, described shortly, the self-interaction term is a constant of the capacitance work function, and does not contribute to the linkage energy.

In principle, one can systematically perform Hill plot analysis with markers for each type of particle in a regulatory network, thereby constructing a complete allosteric map of the protein (see Fig. 11). The problem of optimizing network parameters (interaction energies and particle displacements) from such a dataset in the face of incomplete or noisy data will not be dealt with here, except to say that in rare cases the solution may not be completely constrained, even under ideal conditions. In a unimodal system (responsive to a single external force), with multiple copies of each particle species, it may not be possible to distinguish between homologous and heterologous coupling energies, even when taking advantage of reciprocal relations. In such cases, which will not be considered further here, reciprocal Hill relations between two particle species yield two linkage equations (one for each species) and three unknowns (two self-interaction energies and one allosteric energy).

### Scheme 2 (the HA model)

Having just derived several linkage relations for Scheme 1, we now acknowledge its inadequacy as a model for BK. In reality, the majority of gating charge movement in BK is thought to occur in specialized voltage sensors (J), as indicated by the observation that the Q-V curve is variably steeper and left-shifted compared with the G-V curve (Horrigan and Aldrich, 2002). Adding an additional four J particles, one per subunit, to Scheme 1 generates the nine-particle HA model, which we refer to as Scheme 2 (Fig. 2 A). This is represented by the following linkage diagram (one of six possible diagrams for Scheme 2):

(SCHEME 2)

We see that in upgrading from Scheme 1 to Scheme 2, we have increased the number of allosteric interactions from one (C) to three (C, D, and E). Specifically, pore opening (L) increases the equilibrium constants J and K by factors of C and D, respectively, and J particle activation increases the equilibrium constant K by a factor of E.

The partition function for Scheme 2 is again described by Eq. 12, but the CPFs ZC and ZO, derived from the parenthesized contents in the linkage diagram, are now functions of J and K (Table 4). A single polynomial function f1 describes both CPFs: ZC = f1(J, K, E) and ZO = f1(JD, KC, E). Depending on whether J or K is varied, the CPF ratio ZO/ZC can be interpreted as an expectation value for either D or C. In analogy to Eq. 13, we have:

$ZOZC = f1(JD,KC,E)f1(J,K,E),$
(22a)
$(ZOZC)J(±)= 〈C〉K,$
(22b)
$(ZOZC)K(±)= 〈D〉J,$
(22c)

where it is assumed that E is constant. Eq. 22b corresponds to the sub-scheme in which the voltage-sensor J is maintained (fixed) in an extreme position, allowing K–L linkage to be evaluated, whereas Eq. 22c describes J–L linkage with K fixed. Recognizing from Table 2 that, for BK, JΔ↔VΔ and KΔ↔μΔ, we can apply the respective lever operators to Eq. 15, yielding the following linkage equations:

$ΔμWH[g]v(±) = −4ΔWC,$
(23a)
$Δv(WH[g]+ηL)μ(±)=−4ΔWD.$
(23b)

We should note the somewhat subtle distinction (apart from the superfluous ηL term) between Eq. 16 (Scheme 1) and Eq. 23a (Scheme 2). The latter specifies that voltage be held either at V(−) or V(+), lest the mean activation state of voltage-dependent J particles influence the evaluation of L–K coupling through secondary linkage. We should note that secondary linkage occurs only if J is simultaneously coupled to both L and K particles, i.e., if both ΔWD and ΔWE are nonzero. Otherwise, J could be viewed as merely an accessory component of either the pore or the Ca2+ sensor.

Eqs. 23a and 23b describe how two out of the three allosteric factors in Scheme 2 (C and D) can be derived from Hill transformation of the G-V plot. The third coupling factor, E, is inaccessible to conductance Hill analysis because it contributes equally to ZC and ZO (see Eq. 22a and Table 4). Therefore, E must be determined from the Q-V plot. But before discussing global linkage analysis, we must first finish dissecting the Scheme 2 partition function.

### Higher-order parsing of the Scheme 2 partition function

It is useful to subdivide (parse) the Scheme 2 CPFs ZC and ZO with regard to the resting (R) and activated (A) states of the voltage-sensing J particles, as follows:

$ZC=ZCR+JZCA,$
(24a)
$ZO=ZOR+JDZOA.$
(24b)

The sub-schemes for the possible configurations of voltage-sensing L and J particles are given by the four (Ca2+-dependent) CPFs seen on the right side of Eqs. 24a and 24b, described in Table 5. Note that these second-order CPFs, like their first-order counterparts ZC and ZO, share a single polynomial function, in this case: f2(λK) = 1 + λK, where λ = 1, C, E, or CE. By now, it should be apparent that judiciously chosen ratios of these four CPFs will yield expectation values for C and/or E. For example, Eq. 22b could be rewritten as ZOR/ZORZOA/ZOA = 〈Cμ. Alternatively, ZC and ZO could be parsed with respect to the free (F) and bound (B) states of the K particle, yielding a second set of second-order CPFs of the form g2(λJ) = 1 + λJ, where λ = 1, D, E, or DE (see Table 6). Rephrasing Eq. 22c, we obtain ZOF/ZOFZOB/ZOB = 〈DV. Either parsing route, if connected to a suitable work function, also enables one to measure ΔWE, as demonstrated shortly using a J-based method.

A third and final order of parsing generates CPFs for the eight permutations of (L, J, K) activation states in Scheme 2 (Table 7). In the absence of a fourth regulatory particle added to Scheme 2, these third-order CPFs all equal 1. However, the corresponding LPFs are logarithmically related to the energies of the terminal states, and when entered into a linkage cycle, they provide the correct interaction energies. For example, to determine ΔWE, one could consider the following activation cycle in which J and K are principal and secondary activators, and L is maintained in either the closed or open state:

$ΔWE= ΔKΔJΦL[±]=−kTln[(ZL(±)J(+)K(+)ZL(±)J(−)K(+))/(ZL(±)J(+)K(−)ZL(±)J(−)K(−))].$
(25)

Although not helpful in determining the intermediate points in a linkage plot, the fully parsed states represent the eight “corners” of the five linkage cycles generated by the G-V and Q-V curves in Fig. 2 B. Four of these cycles have already been characterized in the form of Eqs. 23a and 23b. The fifth cycle, shown in the center of Fig. 2 B, measures the interaction energy between charge-carrying and Ca2+-binding elements, and is described next.

### The electrical capacitance energy WC[q]

We recall that the electrical capacitance energy WC[q] is the work function related to the Q-V curve. It is defined as the negative work −Wq of displacing total gating charge per channel (ΔqT) and is equal to the thermodynamic integral:

$WC[q]=−∫0ΔqTVd〈q〉.$
(26)

Integrating Eq. 26 by parts, we obtain:

$WC[q]=∫V(−)V(+)〈q〉dV−ΔqTV(+).$
(27)

Although both terms in Eq. 27 diverge, WC[q] converges to a finite value for well-behaved Q-V plots (i.e., those that approach a maximal value). This is easily seen in the graphical interpretation of Eq. 27 provided by Fig. 3 A. The integral term on the right side of Eq. 27 is equal to A1 + A3 in Fig. 3 A. The product ΔqTV(+) is equal to A2 + A3. The difference between these two areas, −(A2A1), is the negative integral of the Q-V curve reflected across the 〈q〉 = V axis, which precisely describes Eq. 26.

Other useful geometric interpretations of WC[q] make use of the midpoint voltage VM, which was defined in introductory remarks as WqqT = −WC[q]qT. From Eq. 26, we can write the following integral expression for VM:

$VM= ΔqT−1∫0ΔqTVd〈q〉.$
(28)

Chowdhury and Chanda (2012) recently introduced VM to the ion channel literature and named it the “median voltage of charge transfer” in reference to Wyman’s analogous use of the term “median ligand activity” (Wyman, 1964). The word “median” is an apt geometric descriptor, given that a vertical line passing through VM divides a Q-V plot into equal areas, even when normalized to a maximum value of 1.0 (Fig. 3 B). The derivation is as follows: the area of the normalized Q-V in Fig. 3 B that corresponds to the shaded area in Fig. 3 A is A6 + A7A4, which represents the value of WqqT and is therefore equal to VM. In turn, VM is also equal to the area A5 + A6. Equating the two areas, we obtain A7 = A4 + A5, which, as anticipated, divides the shaded regions in Fig. 3 B into equal parts. The fact that VM can be obtained from the unit-normalized Q-V curve is significant, because experimentally it is difficult to measure the “single-channel” Q-V curve (Fig. 3 A), whose maximum value is ΔqT = Qmax/N. Assuming that the value of ΔqT is constant (later relaxed in the setting of variable interaction energies), one would need to measure it only once, using a so-called “charge per channel” experiment (Aggarwal and MacKinnon, 1996; Seoh et al., 1996). The more easily acquired VM (which is independent of N) can be used to monitor the work function WC[q] = −ΔqTVM.

Despite the historical precedent for using the adjective “median” to describe VM—arguments for which include not only the “equal areas” property but also the fact that when V = VM, there are equal populations of channels in the least and most saturated charge states (see Chowdhury and Chanda, 2012)—work by Di Cera and Chen (1993), studying ligand-binding proteins, provides a compelling argument for the use of “mean” instead of “median.” Changing the integrand of Eq. 28 from 〈q〉 to V, one obtains an alternative expression for VM as a function of the normalized electrical capacitance fq = ΔqT−1dq〉/dV (Fig. 3 C):

$VM=∫V(−)V(+)VfqdV.$
(29)

Di Cera and Chen (1993) point out that capacity functions like fq are distributions whose moments—the first of which in an electrical system is VM, and which in a ligand-binding system is μM—are intimately related to coefficients in the partition function.2 Adding the fact that, as we discuss next, there are several practical methods for measuring fq that make it easy to compute VM by using Eq. 29, there are compelling reasons to use “mean” in describing VM, although the historical use of “median” is certainly not wrong for the reasons cited. Here, we will refer to VM as the “mean voltage of activation.”

Of the methods used to acquire fq, probably the least efficient is measuring the slope of the normalized Q-V plot, generally obtained by integrating a series of gating currents obtained by stepping the voltage to different values. More direct methods of measuring fq are possible, including computing admittance for noise-driven gating currents at different holding potentials (Fernández et al., 1982), or measuring the gating current in response to a slow voltage ramp (Sigg and Bezanilla, 1997). The voltage-ramp method, which has the advantage of generating the entire voltage dependence of fq in one sweep, is later demonstrated here through simulation. We note that fq is proportional to the equilibrium variance of gating charge displacement and is therefore positive for all potentials (obviously a necessary property for a distribution function).

A final relevant property of VM in systems that are both voltage and ligand sensitive is its rate of change with respect to ligand activity. To derive this, we differentiate Eq. 28 with respect to μ to obtain:

$ΔqTdVMdμ=∫0ΔqT(∂V∂μ)qd〈q〉=−∫0ΔqT(∂〈n〉∂〈q〉)μd〈q〉,$
(30)

where the Maxwell relation (∂V/∂μ)q = −(∂〈n〉/∂〈q〉)μ is used. After evaluating the integral on the right, we have:

$ΔqTdVMdμ=− Δv〈n〉,$
(31)

where VΔ〈n〉 is the change in the number of bound ligand 〈n〉 in response to a maximal change in voltage. Measuring the maximum value of dVM/ in BK channels, and knowing ΔqT, one can obtain a lower bound (= VΔ〈nmax) for the total number of calcium-binding sites ΔnT as an alternative to measuring the Hill coefficient.

### A statistical mechanical formulation of WC[q]

Recalling that 〈q〉 = −(∂Φ/∂V)μ, and invoking the “bridge” equation (Eq. 4), we derive from Eq. 27 the following:

$WC[q]= kTln(ZV(+)ZV(−))−ΔqTV(+).$
(32)

The LPFs ZV(−) and ZV(+) in Eq. 32 correspond to sub-schemes of Scheme 2 derived from extreme potentials V(−) and V(+). ZV(−) describes J and L particles as fixed in their resting positions, whereas K particles are allowed to bind Ca2+, and ZV(+) is analogous except with J and L fixed in their activated positions. The LPFs contain CPFs ZCR and ZOA, as follows: ZV(−) = ZCR4 and ZV(+) = LV(+)(JV(+)DZOA)4. Recalling that ηL = −kTlnL and ηJ = −kTlnJ, we insert the above expressions for ZV(−) and ZV(+) into Eq. 32, yielding:

$WC[q]= 4kTln(ZOAZCR)−ηL,V=0−4ηJ,V=0−4ΔWD,$
(33)

where we have eliminated the diverging second term in Eq. 32 by recognizing that, because ΔqT = ΔqL + 4ΔqJ, the electrical components of ηL,V(+) and ηJ,V(+) exactly cancel ΔqTV(+), leaving only the charge-neutral potentials ηL,V=0 = ΔGL − ΔnLμ and ηJ,V=0 = ΔGJ − ΔnJμ. In unusual cases where obligatory calcium binding is associated with the activation of L and/or J particles, one would observe sloping asymptotes (m = ΔnL + 4ΔnJ) when plotting WC[q] versus μ, but for standard models of BK, the global linkage plot has horizontal asymptotes.

The sigmoidal component of WC[q] versus μ in Scheme 2 possesses the vertical span:

$ΔμWC[q]=−ΔqTΔμVM=−4(ΔWC+ΔWE).$
(34)

The right side of Eq. 34 follows immediately from Eq. 33 by substituting the expression (1 + CEK)/(1 + K) for ZOA/ZCR, and recognizing that this is equal to the expectation value 〈CEK. The outcome of Eq. 34 solves the earlier problem of not being able to measure E solely through conductance Hill analysis, because now the value of 4ΔWE can be obtained by subtracting Eq. 34 from Eq. 23a.

### Limited global analysis

It is reasonable to enquire, given the outcome of Eq. 34, whether a limited version of global analysis can be used to measure ΔWE independently of ΔWC. The answer is a qualified yes, provided one can maintain the pore in either the closed (L(−) = C) or open (L(+) = O) state, thereby eliminating its influence on J–K linkage. With the pore locked in this manner, Eq. 34 can be rewritten:

$ΔqJΔμVM[C or O] = −ΔWE,$
(35)

where VM(C) and VM(O) are mean activation potentials derived from the limiting curves 〈qC〉 = kT(∂lnZL(−)/∂V)μ and 〈qO〉 = kT(∂lnZL(+)/∂V)μ. Using Eq. 12, these can be expanded to:

$〈qC〉= 4kT(∂lnZC∂V)μ,$
(36a)
$〈qO〉= ΔqL+4kT(∂lnZO∂V)μ.$
(36b)

In BK channels, it is possible to resolve 〈qC〉 and 〈qO〉 through kinetic means by taking advantage of the fact that pore opening is ∼100 times slower than voltage-sensor activation (Horrigan and Aldrich, 2002).

Purely thermodynamic methods can also be used to determine 〈qC〉 and 〈qO〉 from measurements of G-V and Q-V plots, for example, through the relations 〈qC〉 = 〈q〉 + kTdln(GmaxG)/dV and 〈qO〉 = 〈q〉 + kTdlnG/dV. The expression kTdlnG/dV is known as the mean activation charge displacement 〈qa〉, first used in reference to the Shaker K+ channel (Sigg and Bezanilla, 1997), where it appears to follow the course of an “upside-down” Q-V curve. In BK, measurements of 〈qa〉 (Horrigan and Aldrich, 2002) are consistent with an allosterically regulated pore particle, as its value at very negative voltages (the so-called “limiting slope”) decreases to a fairly small value (∼0.3 eo), compatible with that of ΔqL. In contrast, the limiting slope observed in Shaker is substantially larger (12–13 eo; see Islas and Sigworth, 1999)—equal in value to ΔqT measured from “charge per channel” experiments (Aggarwal and MacKinnon, 1996; Seoh et al., 1996)—consistent with a single open state at the end of a voltage-dependent activation sequence.

A different but related method uses the slope of the conductance Hill to determine 〈qC〉 and 〈qO〉. Starting from the identity WH[g] = kTln(ZL(+)/ZL(−)), we derive dWH[g]/dV = 〈qO〉 − 〈qC〉 (see also Conti, 1986), and recognizing that 〈q〉 = (1 − PO)〈qC〉 + POqO〉, we obtain 〈qC〉 = 〈q〉 − (G/Gmax)dWH[g]/dV and 〈qO〉 = 〈q〉 + (1 − G/Gmax)dWH[g]/dV. With both of these methods, as with the earlier use of Eqs. 23b and 34, a combination of global (Q-V) and local (G-V) techniques is necessary to measure ΔWE for the simple reason that the J–K interaction does not directly involve the pore.

### Model-independent application of linkage analysis

Linkage plots yield model-independent information about channel gating, starting with whether the channel is part of an “accessible” allosteric network. We have just discussed how 〈qa〉 experiments support the notion of obligatory voltage-sensor activation before the opening transition in Kv channels like Shaker. This would constitute a “tight” allosteric network in which the coupling energy (which could be negative, as from steric interference) between voltage sensors and the pore is so large as to be immeasurable using electrophysiological technique. Models of Shaker proposed to date do not predict conductance Hill plots that resemble Fig. 1; rather, one expects dramatically different slopes for negative (mV(−) = ΔqT) and positive (mV(+) = ΔqL) asymptotes.

In nonobligatory or accessible allosterism, work functions related to identifiable structural components can be used to construct a network map consisting of particle displacements and interaction energies, with a view of understanding structure–function relationships from an energetic standpoint. Turning again to the BK channel, we would like to develop linkage relations that are not confined to specific models like Schemes 1 and 2, but relate to more general schemes that can be constrained through experiment but are imprecise with respect to internal processes. The partition function Z is a useful tool in implementing such “fuzzy” modeling. We have already demonstrated its utility in describing simple models. The exact partition function of a protein cannot be known exactly, but one can take advantage of the fact that Z is essentially a sum of probabilities and can be parsed into LPFs representing experimentally distinguishable limiting states. In the case of the BK channel, the usual modeling constraints are as follows, supported by numerous studies (Cox et al., 1997; Rothberg and Magleby, 2000; Horrigan and Aldrich, 2002; Xia et al., 2002):

#### (1) Voltage-sensing and calcium-sensing components of the channel are distinct.

Specifically, the voltage sensor and pore are both voltage sensitive (ΔqL, ΔqJ ≠ 0), whereas the calcium-binding domain demonstrates no voltage sensitivity (ΔqK = 0), a distinction made in Table 2. This is a reasonable assumption. Both the pore and voltage sensor are located in the membrane, where charged residues are sensitive to the membrane electric field. The known calcium-binding elements in the gating ring are located in the cytoplasm, displaced internally from the membrane plane, and therefore calcium binding is not expected to traverse a significant portion of the membrane potential, minimizing its contribution to “gating charge.”

#### (2) Pore conductance is binary, transitioning stochastically between closed (C) and open (O) states.

This implies that the pore can be modeled as a two-state process. The experimental proof in BK, as in other channels, is found in tracings of single-channel ionic currents, which undergo very rapid transitions between nonconducting and conducting states. Very detailed analysis of single-channel ionic currents cast doubt on the conjecture that pore conduction is a purely binary process, as subconductance states in BK have been observed under various conditions (Stockbridge et al., 1991; Ferguson et al., 1993; Mistry and Garland, 1998; De Wet et al., 2006). The consequences of a nonbinary pore will be addressed later. For now, we accept the binary pore as a given.

#### (3) Channel subunits do not interact except through the pore.

BK channels are composed of four identical subunits centrally connected to the pore domain. It is usually assumed that contact between subunits is minimal except where they interact with the pore. This is clearly plausible for the four voltage-sensor components, which are separated from each other in space (Wang and Sigworth, 2009). It is not as clear whether intersubunit interactions occur in the gating ring, where there is extensive contact between calcium-binding domains from different subunits (Yuan et al., 2010). Nevertheless, experiments suggest that the gating ring is dominated by intrasubunit, not intersubunit, interactions (Qian et al., 2006). This assumption simplifies Z by allowing the subunit CPF to be raised to the fourth power. On the other hand, the existence of intersubunit interactions can be easily dealt with by expanding Z with respect to powers of the equilibrium constant of the particle of interest, as demonstrated earlier for Scheme 1 and again later for the multi-subunit pore (Scheme 5).

#### (4) Allosteric energies are constant, independent of voltage or calcium concentration.

This last point is the weakest in terms of experimental evidence supporting it. There is no a priori reason why allosteric interactions should be independent of environmental influence. It has been postulated, for example, that interactions between voltage-sensing particles in the voltage sensor might also be voltage sensitive (Pantazis et al., 2010), and interactions between the gating ring and the voltage sensor appear to be mediated by internal Mg2+ (Shi et al., 2002; Yang et al., 2008). The existence of environmentally sensitive allosteric factors distorts linkage relations, as discussed later.

These constraints on BK function can be translated into a general gating scheme (Scheme 3), of which the HA model (Scheme 2) is the simplest representation. Scheme 3 (illustrated by the cartoon in Fig. 4) possesses a binary pore, but the activation of other modality-specific domains may be more complex, characterized by composite equilibrium constants (J, K) and internal allosteric factors (F, G), as indicated by bold type. The linkage diagram of Scheme 3 in the parsing order LJ[G]K[F] is given by:

(SCHEME 3)

The resemblance to Scheme 2 is obvious (in particular, Eq. 12 still applies), with the difference lying in the ambiguous nature of the number and configuration of particles within J[G] and K[F]. Nevertheless, it is possible to mathematically characterize Z and derive linkage relationships. The details of these derivations will not be presented here but can be found in the supplemental text. Not surprisingly, the linkage equations arising from such a process are similar to those already derived for Schemes 1 and 2. They are given by:

$ΔμWH[g]v(±) = −4ΔWC,$
(37a)
$Δv(WH[g]+ηL)μ(±)=−4ΔWD,$
(37b)
$ΔμWC[q]= −ΔqTΔμVM = −4(ΔWC+ΔWE),$
(37c)
$ΔμWC[q]L(±)= ΔqJΔμVM[C or O] = −4ΔWE,$
(37d)

where we are reminded that the energies and displacements of composite variables are the sums of their constituents. One important point to recognize is that Eqs. 37ad, which are derived from work functions related to the G-V and Q-V curves, do not specify the value of the internal allosteric factors F and G. To probe the inner operation of either the gating ring or voltage sensor using the Hill analysis, one requires specific probes for individual particles within these domains.

This concludes the theoretical section of the paper. A summary of linkage analysis formulas used for the generalized Scheme 3 and its derivatives with regard to G-V and Q-V curves is given in Table 8. To investigate the challenges in constructing and analyzing linkage plots, simulations of gating and ionic currents from a “complex” 17-particle kinetic model (Scheme 4) were performed under patch-clamp–like conditions, as described previously in Materials and methods. The thermodynamics of Scheme 4 and the results of these simulations are described next.

### Thermodynamics of the 17-particle model (Scheme 4)

The test model (Scheme 4) used in kinetic simulations is shown in Fig. 5. The scheme consists of a central pore (L) surrounded by four identical regulatory subunits, with each subunit containing two voltage-sensing particles (J1 and J2) and two calcium-sensitive particles (K1 and K2), for a total of 17 interacting gating particles. The particles L, J2, K1, and K2 form an allosteric loop involving the core interactions C, D, and E, and an internal interaction F coupling K1 and K2. The J1 particle does not participate in the loop but influences its neighboring J2 particle through an internal factor G. Although this scheme has not been experimentally validated, the particle arrangements and values of model parameters represent an amalgam of interactions previously proposed for BK (see references in Table 9). The linkage diagram for Scheme 4 parsed in the order L→J2→(J1, K2)→K1 is given by:

(SCHEME 4)

from which we derive, as usual, the first-order equation Z = ZC4 + LZO4. Expressions for the CPFs ZC and ZO are as follows:

$ZC= (1+J1)ZCR+J2(1+J1G)ZCA,$
(38a)
$ZO=(1+J1)ZOR+J2D(1+J1G)ZOA.$
(38b)

The second-order CPFs ZCR, ZCA, ZOR, and ZOA are functions of K1 and K2, given by ZCR = f2(K1, K2) = (1 + K1) + K2(1 + K1F), ZCA = f2(K1E, K2), ZOR = f2(K1, K2C), and ZOA = f2(K1E, K2C). The mean single-channel gating charge displacement 〈q〉 ranges from 0 to ΔqT = ΔqL + 4(ΔqJ1 + ΔqJ2) and is calculated from 〈q〉 = 〈l〉ΔqL + 〈j1〉ΔqJ1 + 〈j2〉ΔqJ2, where the voltage-dependent particle activation curves are derived as follows:

$〈l〉≡∂lnZ∂lnL=LZO4Z,$
(39a)
$〈j1〉≡∂lnZ∂lnJ1=(4J1Z)[(ZCR+J2GZCA)ZC3+L(ZOR+J2DGZOA)ZO3],$
(39b)
$〈j2〉≡∂lnZ∂lnJ2=(4J2Z)[(1+J1G)ZCAZC3+LD(1+J1G)ZOAZO3].$
(39c)

The average single-channel conductance is 〈g〉 = 〈lgL.

### Linkage analysis of simulated data using Scheme 4

Numerically computed equilibrium curves derived from Scheme 4 are shown in Figs. 6 (WH[g] “Hill” analysis) and 7 (WC[q] “capacitance” analysis) for a range of [Ca2+]i values. Because Scheme 4 is a special case of the generalized Scheme 3, Eqs. 37ad were used to numerically compute work functions. Monte Carlo–simulated currents were analyzed according to the procedure described in Materials and methods. The simulation results are shown in Figs. 810, with averaged parameters from 10 experiments listed in Table 9. Despite Nyquist and channel fluctuation noise, parameter estimates were fairly accurate, with a maximum error of 7.2% in the global parameter ΔWC + ΔWE, and three out of five averaged outcomes found to be statistically indistinguishable from their true values.

The “charge-per-channel” experiments (Fig. 8), in which N was estimated using ionic mean–variance analysis and Qmax was obtained by integrating gating currents, slightly underestimated the true value of ΔqT = 2.62 eo by −2.4%, a difference that is close to being statistically significant (P = 0.055).

The major challenge in generating WH[g] plots from simulated data (Fig. 9) was determining the “floor” and “ceiling” of the G-V curve. These limits depended on the randomized “patch-specific” values of N, Veq, and gleak. From the definition WH[g] = kTlnGkTln(GmaxG), we recognize that the kTlnG defines the V(−) asymptote of WH[g] and is logarithmically sensitive to the floor, whereas kTln(GmaxG) defines the V(+) asymptote and is logarithmically sensitive to the ceiling. Although Nyquist noise may have played a small role in the uncertainty of the G-V near the asymptotes, a far greater source of variation was open–closed fluctuations. The signal to noise (S/N) ratio for the floor variable G from N independently gating channels is derived from the mean G = NgLPO and variance VarG = NgLPO(1 − PO) through the relation S/N = G/VarG1/2, which works out to (NG/(GmaxG))1/2. For S/N > 1, many of the lower asymptotic points kTlnG were undefined, because fluctuations from the mean resulted in zero or even negative conductance (recall that the baseline or “floor” of the conductance was not known at the time of analysis and needed to be determined). This condition occurs when PO < 1/N. Therefore, with a cumulative N of 105 (= 102 traces × 103 channels), noise effects became evident for PO ≤ 10−5, as seen in Fig. 9. It is remarkable that useful data could be derived for PO smaller than the predicted limit, down to PO = 10−8. The corresponding S/N ratio for the ceiling variable GmaxG is (N(GmaxG)/G)1/2. Thus, for S/N ≈ 1, and cumulative n = 105, noise is increased when 1 − PO = 1/N, corresponding to the upper limit of the plot in Fig. 9. Because the V(−) asymptote at high calcium concentration (μ(+)) occurred at intermediate values of PO, it was least affected by floor and ceiling effects and was used to fit the slope ΔqL. The average fitted value of ΔqL differed only slightly (+1.6%) from the given value of 0.3 eo. To determine interaction energies ΔWC and ΔWD, the relative heights of the other three “noisy” asymptotes were estimated through a process of alternatively fitting and manually adjusting the patch-specific variables (N, Veq, and gleak) so that they visually paralleled the (V(−), μ(+)) asymptote. The value of ΔWC (averaged from measurements taken at V(−) and V(+)) was 1.3% less than the given value of −53 meV. Measuring ΔWD was slightly less precise, with a −3.0% underestimation of the given value of −81 mV (P = 0.009), decreasing to only −0.7% if the noisier “zero” Ca2+ plot was excluded from the analysis.

In WC[q] experiments (Fig. 10), the average value of the outcome parameter ΔWC + ΔWE differed significantly from the given value of −75 meV by −7.2% (P = 0.044). Nyquist noise was a major factor, with simulated gating currents barely exceeding 15 fA compared with the root-mean-square noise level of 29 fA at 1 kHz. After signal averaging 100 traces and post-filtering to 0.1 kHz, a noise level of 0.9 fA was reached, which improved signal recognition but did not alleviate the uncertainty in assigning proper baselines for integration, which had to be performed manually by eye. In the final analysis, however, a significant source of error in WC[q] analysis was not random but systematic, caused by a “lag” in the gating current from failing to maintain equilibrium near the steep portion of the Q-V curve, where the channel responds mostly slowly. This effect was more pronounced in high than in low [Ca2+]i. It was quantified by simulating the model in the absence of Nyquist noise, and the VΔVM of the now smooth Q-V curves (shown in Fig. 10 as thick lines) differed from the true value by −6.9%. Thus, the lag effect, when added to the small inaccuracies in ΔqT measurements, accounted for most of the error in ΔWC + ΔWE.

### Hill analysis performed on markers of voltage- and Ca2+-sensor activation

Hill analysis is not limited to G-V plots, provided that markers of activation can be obtained for particles other than the pore. Energy Hill plots for each of the four subparticles (J1, J2, K1, and K2) in Scheme 4 are shown in Fig. 11. The work functions for the Ca2+-sensitive K particles are plotted on μ and log[Ca2+]i axes, replicating traditional ligand-binding curves. The displacements and interaction energies derived from these linkage plots are easily matched to the highlighted parameters shown in the cartoons (Fig. 11). The rise of the sigmoidal components with respect to the slope of the asymptotes varied among particles, with J1 providing the best resolution because of its small intrinsic charge and relatively large interaction energy with J2. From an experimental standpoint, the Hill plot for J2 would be considered the most problematic among the voltage-sensing particles, because limiting behavior was not reached even for very large values of V, leading to the risk of measuring incorrect slopes and interaction energies from “false” asymptotes. The Hill plots for the K particles (Fig. 11, C and D) are notable in demonstrating net negative displacements in the sigmoidal component stemming from negative cooperativity between K1 and K2WF > 0).

### Adding a tetrameric pore: The 20-particle model (Scheme 5)

We now reexamine the earlier constraints (1–4) on BK gating by relaxing them in Scheme 4. Constraint 2 posited that pore activation is binary, based on single-channel recordings. However, this is inconsistent with the tetrameric structure of the pore and does little to explain the presence of observed short-lived subconductance events. A more compelling hypothesis has been put forward (Chapman et al., 1997; Zheng and Sigworth, 1997), which postulates that the four P subunits, through strong homologous interactions, activate in near concerted fashion, allowing only brief visits to intermediate conducting states. A modification of Scheme 4 that includes such a tetrameric pore is shown in Fig. 12 A (Scheme 5).

As was done previously for K particles in Scheme 1, multiplicative constants B1, B2, B3, and B4, in this case representing nearest-neighbor P-particle interactions (but excluding trans-interactions), were used to modify the equilibrium constants P of successively activating pore particles. The intrinsic P particle transition energy ΔGP = ΔGL/4 was additionally offset by the amount −ΔWB = −(ΔWB1 + ΔWB2 + ΔWB3 + ΔWB4)/4 to maintain the original free energy change ΔGL between closed and open states. Under these conditions, of the six configurational pore states (which include the cis- and trans-configurations of the p = 2 state), only the closed and open states experience a match between the number of interactions and activated particles (Fig. 12 B). Intermediate (and presumably subconducting) states are characterized by an imbalance of these numbers, leading to increased state energies (on the order of −ΔWB and −2ΔWB; Fig. 12 B) and brief occupancy times. These attributes of the interacting tetrameric pore are found in the Scheme 5 partition function:

$Z = ZC4+4BZC3ZOP+(2+4B1B2)ZC2ZC2P2+4B1B2B3ZCZO3P3+ZO4P4,$
(40)

where ZC and ZO remain unchanged from Scheme 4 (Eq. 38). It is readily apparent that for large values of B = exp(−ΔWB/kT), the partition function reverts back to Scheme 4 with its binary pore equilibrium constant L = P4. Although we have assumed that P–P interactions favor opening transitions, we could have easily constructed a scheme where the activation pathway runs in reverse, with positive interactions favoring closing transitions, as long as the individual offset energies are adjusted to maintain the desired equilibrium between closed and open states.

The effect of a strongly cooperative multi-subunit pore on work functions is shown in Fig. 12 C for ΔWB = −250 meV. Assuming equal increments of interaction energies and conductance increases in the pore activation sequence p = {1…4}, model parameters were determined as follows: Bp = B1/4, 〈g〉 = 〈pgL/4, 〈p〉 = ∂lnZ/∂lnP, and ηP = ηL/4. The results of global capacitance (WC[q]) analysis were practically unchanged from Scheme 4 because interactions between voltage- and Ca2+-sensitive particles are unaffected by B; therefore, WC[q] plots are not shown. However, dramatic decreases in the limiting slopes of WH[g] are demonstrated at very extreme voltages (Fig. 12 D), whereas for a more attainable range (−100 to 300 mV), the behavior was nearly in line with the usual (Scheme 4) predictions of a binary pore. The extreme V asymptotes of the Hill function WH[g] were equal to:

$WH[g]V(−)=ΔWB+kTln(ZORZCR)−ηP,V(−),$
(41a)
$WH[g]V(+)=−ΔWB−ΔWD+kTln(ZOAZCA)−ηP,V(+).$
(41b)

The rise VΔ(WH[g]ηP)μ(±) in asymptotes taken from the difference of Eqs. 41a and 41b equals −(ΔWD + 2ΔWB). The value of −2ΔWB for the self-interacting term is consistent with the earlier statement that the self-interaction term is equal to the difference in cooperative energies of the first and last transitions. It is appreciated that the difficulty in measuring ΔWB in native BK channels would be considerable, considering the large potentials needed to reach “true” asymptotic behavior and the need to resolve subconductance levels.

In Table 8, equations are provided for the general case, which is characterized by arbitrary values for the self-allosteric factors Bp and subconductance levels gp = xpgL. The loss of uniform incremental values in these parameters influences the measurement of self-interaction energy, as demonstrated by the new linkage relation:

$Δv(WH[g]+ηP)μ(±)=−(ΔWD+ΔWB3+ΔWB4)+kTln[16x1(1−x3)].$
(42)

The second term on the right of Eq. 42 is an artifact of the nonuniform spacing of subconductance levels and vanishes if x1 = 0.25 and x3 = 0.75. In an apparent contradiction with the earlier assertion that the self-interaction energy is a function of the first and last steps of activation, Eq. 42 appears to depend on interaction energies from the last two steps. The contradiction resolves when one considers that, according to the partition function in Eq. 40, the first and last steps contribute the following energies: kTln(1/B) and kTln(B3/B1B2), the difference of which is kTln(B3B4) = −(ΔWB3 + ΔWB4). These expressions stem from the constraints imposed on self-interaction energies to achieve a “balanced” equilibrium between closed and open states.

### Dual-modality particles

Constraint 1 for BK gating required all gating particles to be unimodal, i.e., each responsive to only one external force. However, it is not hard to imagine circumstances in which this would be false. For example, a Ca2+-binding site might be displaced from the internal solution by a fraction d of the membrane potential, leading to a gating charge displacement of 2d eo. Another mechanism for dual sensitivity is the obligatory binding of Ca2+ during pore opening. As unlikely as this seems, Piskorowski and Aldrich (2002) have reported Ca2+ sensitivity in a BK channel whose gating ring had been removed, suggesting that removal of this bulky domain opened up alternative binding sites, perhaps even on the pore itself. A modification of Scheme 4 in which the pore partially binds Ca2+ upon opening was performed by specifying the particle function μL = ΔGL + −ΔqLV − ΔnLμ, where ΔnL was chosen to be 0.4.3 This led to sloping asymptotes in the WC[q] versus μ plot (Fig. 13 A). The energy Hill plot WH[g] demonstrated asymptotic slopes for both V and μ axes (Fig. 13 B).

### Modality-sensitive allosteric factors

Constraint 4 implied that allosteric energies are insensitive to external forces. There is reason to believe that this might not be the case in BK. For example, Mg2+ appears to play a role in the interaction between gating ring and voltage sensor (Yang et al., 2008), raising the possibility that a small gating charge associated with Mg2+ displacement could affect the interaction energy ΔWE. Another example is a small gating charge that appears to mediate the interaction (described by the internal coupling factor G in Scheme 4) between putative voltage-sensing particles, as inferred from fluorescence experiments (Pantazis et al., 2010). To incorporate voltage sensitivity into allosteric factors, an electrical term can be added to the interaction energy, for example ΔWG = ΔGG − ΔqGV. A voltage-sensitive G = exp(−ΔGG/kT) would contribute an amount (∂lnZ/∂lnGqG of gating charge to the Q-V curve. The effects on the Q-V and Hill plots from a small charge (ΔqG = 0.2 eo) are shown in Fig. 14 (A and B). Because interactions between voltage- and Ca2+-sensitive components of the channel do not depend on internal energies within the voltage sensor, global analysis (Fig. 14 A) was unaffected except that ΔqT grew larger by the amount 4ΔqG. On the other hand, the V asymptotes in the local Hill plot for J2 diverged because of the addition of ΔqG to the slope of the upper asymptote (Fig. 14 B). The interaction energy −(ΔWG + ΔWD) could still theoretically be obtained by evaluating the height differences between asymptotes evaluated at V = 0.

In a second modification to Scheme 4, a charge ΔqE = 0.2 eo was assigned to the allosteric factor E, which, unlike G, couples particles of different modalities (J2 and K1). Here, a more dramatic change to the Q-V curve occurred (Fig. 13 C). At very low [Ca2+]i, where E ≈ 0, the channel behaved as in the original Scheme 4 with a well-defined V(−), despite the theoretical value of VM being shifted immeasurably far to the right on the V axis. With increasing [Ca2+]i, The E charge component (∂lnZ/∂lnEqE entered the Q-V plot from the right side, with VM eventually reaching a well-defined limiting value specified by V(+). The global energy of interaction −4(ΔWC + ΔWE) can be deduced from the area A2A1 shown in Fig. 14 C. From an experimental standpoint, obtaining the correct result would require that one recognize that ΔqT, the total charge movement per channel, has become a function of [Ca2+]i as a result of linkage between J and K. Complex behavior was also demonstrated in the WH[j2] plot (Fig. 14 D), which contained numerous inflection points, generating “false” and ultimately diverging asymptotes across a wide range in voltage.

### Intersubunit interactions

Finally, we briefly address constraint 3, which stated that subunit interactions between regulatory domains are forbidden except when acting through the pore. In light of generally extensive interactions between cytoplasmic domains found in many K+ channels, known to be important for channel assembly, as reviewed by Schwappach (2008), this may seem unrealistic. It is interesting to note that voltage-sensor domains in Shaker appear to be anchored to neighboring subunits in an effort to increase the mechanical advantage for gating the pore (Long et al., 2005), an arrangement that falls outside the models considered here. As demonstrated in a variation of Scheme 1 and in Scheme 5, the partition function Z for a channel possessing interacting subunits must be expanded in powers of the equilibrium constants of interacting particles. This presents no serious difficulties apart from increasing the complexity of Z compared with Eq. 12. As with independently acting regulatory domains, the primary effect of positively self-interacting or intersubunit-interacting protomers in K+ channels is to strengthen allosteric coupling with the pore, thus steepening the G-V curve and shifting it to the left.

Ion channels achieve specialization by acquiring sensory units that activate in response to external forces (voltage, chemical potential, membrane tension, etc.). The ensuing displacements (charge, ligand binding, strain, etc.) may influence pore opening directly or through allosteric means. A thermodynamic description of a channel’s gating network would include, in addition to the spatial arrangement of regulatory particles, the values of displacements and interaction energies. The dual-regulated BK channel serves as an ideal system in which to study such energy/displacement networks, as it has a large configurational space that appears to be fully accessible to electrophysiological experimentation. Interactions between the pore and regulatory domains in BK remain incompletely understood at a molecular level. The eventual elucidation of a mechanistic description of gating will likely be achieved through careful probing of mutual interactions between specific residues using site-directed mutagenesis or other site-specific techniques. The interaction energies and displacements obtained from linkage analysis are useful metrics by which such experimental interventions can be interpreted. In studying ion channels, the principal sources of linkage information are the Q-V and G-V plots, representing normalized equilibrium curves for the gating charge q and the conductance g, respectively.

In this paper, equilibrium curves derived from gating models of BK were analyzed using two work functions: the electrical capacitance energy WC[q], a global parameter derived from the Q-V plot; and the Hill energy WH, which is a local parameter applied to specific markers of activation, for example the G-V curve in the case of the conductance Hill energy WH[g]. It has been emphasized that linkage analysis is model independent in the sense that only general assumptions regarding channel allosterism are needed to interpret linkage plots. In practice, should the linkage functions behave unpredictably—for example, if the limiting asymptotes of a work function diverge—then the underlying assumption regarding channel gating would need to be reexamined. Otherwise, linkage plots are easy to interpret, with channel-specific displacements and interaction energies equal to the slopes and height differences in V and μ asymptotes. However, noisy data and other artifacts have a corrupting influence, and the reverse problem of constructing the best possible model given a set of possibly unreliable measurements and other constraints should be considered. Although such considerations are outside the scope of this paper, Monte Carlo simulations performed on Scheme 4 yielded some insight into the quality of the data that can be expected from patch-clamp experiments and pointed to some potential obstacles. These are discussed next.

### Monte Carlo simulation of ramp experiments in the patch

The patch-clamp technique, in which a patch of membrane containing one or many channels is electrically isolated using a polished glass pipette to establish a gigaohm seal, is the current method of choice in measuring BK channels because of its excellent performance characteristics and experimental access to the internal solution. The use of slow voltage ramps in the Monte Carlo simulations performed here enabled “quasi-equilibrium” curves to be generated from single sweeps. Patch-clamp conditions were simulated by incorporating Nyquist noise, leak, and linear capacitance currents and randomly assigned values to patch-specific variables to mimic conditions encountered during analysis of real gating and ionic currents.

### Sources of error in electrical capacitance energy measurements

The major source of error in WC[q] measured in voltage-ramp simulations was systematic, caused by a Ca2+-dependent “lag” in gating charge capacitance near the steep portion of the Q-V curve. In principle, this error could be cancelled by reversing the direction of the ramp and averaging the VM values for upsloping and downsloping ramp protocols. Nyquist noise also played a significant role, requiring signal averaging and post-filtering of the simulated gating current. A faster ramp speed would have improved S/N but also increased systematic error. On the other hand, Horrigan and Aldrich (2002) achieved very good S/N in gating capacitance measurements through admittance analysis using a sinusoid stimulus superimposed on a 1-s voltage ramp. Their choice of excitation frequency = 0.868 kHz selectively activated the voltage sensor; however, a white-noise excitation source could theoretically capture the full frequency response of the gating current (Fernández et al., 1982). The definitive solution to noisy data is to increase N with a greater channel density or by measuring larger areas of membrane, for example, through use of the Vaseline gap cut-open oocyte method (Taglialatela et al., 1992), which can record >106 channels and also allows internal perfusion.

When performing global analysis from Q-V curves, care should be taken to rule out sloping asymptotes (observable as VM “creep”) and an environmental dependence in the value of Qmax, as these complicate the analysis (see Figs. 13 and 14). Another potential complication is from nonlinear capacity effects. The ramp elicits a linear membrane capacitance current (= mCpatch) that was assumed to be constant. In reality, there may be electrostriction effects, or an even more insidious process: a small but broad “fast” capacitance linked to individual U-shaped–state potentials in a channel’s conformational energy landscape. In the Shaker channel, the fast capacitance decreases by an amount of ∼5 × 10−4 eo/mV per subunit as the channel transitions from the most closed state to the open state (Sigg et al., 2003). This amounts to only 0.4% of the peak capacitance in Shaker and ∼5% of the peak capacitance predicted by Scheme 4 using parameters in Table 9. It is unknown whether a similar fast capacitance is measurable in BK. If there were a substantial state dependence in the electrical capacitance, say in the form of a difference ΔCJ = ∂ΔqJ/∂V between the resting and activated states of the voltage-sensing J particle, one would need to modify the particle potential of J, which by integrating ∂ΔqJ (assuming ΔCJ is voltage independent) would obtain a second-order voltage term: ηJ = ΔGJ − ΔqoJV − ΔCJ(VVo)V. This could theoretically have an effect on baselines and limits of integration when analyzing ramp-generated Cq to measure VM. Later, we’ll see how a strong state dependence in heat capacity can dramatically influence linkage relations in temperature-activated TRP channels (Clapham and Miller, 2011).

Measuring WC[q] requires knowing the value of total gating charge, the so called “charge-per-channel” ΔqT = Qmax/N. Although a model-dependent estimate of ΔqT = 2.62 eo has been determined for BK (Horrigan and Aldrich, 2002), an independent determination of ΔqT using charge-per-channel methods has yet to be published. Such experiments (see Fig. 8) require that one measure gating currents and N in the same preparation (Aggarwal and MacKinnon, 1996; Seoh et al., 1996).

### Sources of error in conductance Hill energy measurements

The critical step in WH[g] analysis of simulated ionic current data was establishing the “floor” and “ceiling” baselines for the G-V plot, as Hill asymptotes are extremely sensitive to these values. The simulations demonstrated that reliable measurements could be obtained somewhat beyond the theoretical limits of PO < 1/N for the floor asymptote, and 1 − PO = 1/N for the ceiling asymptote.

Several problems predictably arise when considering experimental limitations. One is that the size of BK currents—in the simulations they exceeded 12 nA for a 50-pS pore conductance—would create serious signal distortion caused by series resistance artifact and transient accumulation of K+ ions in the semi-confined space of the outer pore vestibule. A method to reduce current size without altering channel gating would be necessary to accurately measure the large PO (“ceiling”) asymptote of WH[g].

A second experimental issue is the large dynamic range of open probabilities necessary to generate both negative and positive asymptotes in the Hill plot. Scheme 4 simulations required measuring open probabilities as low as 10−7 (see Fig. 9), requiring at the very least a 24-bit A/D converter (107 ≈ 224). An elegant solution to this problem would be to construct a “Hill” circuit, which generates analogue signals for G and Gmax − G by adjusting the driving force VVeq and baseline Ibase of the ionic current, and then applies a logarithmic converter before digitizing. Using a repeating ramp protocol, manual or semi-automated adjustment of inputs Gmax, Veq, and Ibase could be performed to achieve an optimal shape of the Hill curve in real time.

A third problem is experimental noise, which in ionic current simulations primarily consisted of channel-related gating fluctuations rather than Nyquist noise. As discussed earlier, the estimated cumulative number of channels N needed to generate a good quality Hill plot is the inverse of the larger value of PO and 1 − PO. To use markers of activation other than the G-V curve (for example, from fluorescent labels), one would need to consider the S/N characteristics of the signal. Ionic currents related to voltage-sensor activation that could theoretically be used as “markers” in specially engineered ion channels are the “proton pore” current (Starace and Bezanilla, 2004) and “omega” current (Tombola et al., 2005).

### Relation to the HA model

In their comprehensive study of BK gating, Horrigan and Aldrich (2002) used a mixture of kinetic and thermodynamic methods to constrain values of gating parameters in the context of their nine-particle model. The HA model is arguably the first complete description of BK gating and provides the basis for the schemes considered here. Most of the values used in Scheme 4 were derived from the results of this study (Table 9). Many of their methods can also be framed within the broader context of linkage analysis.

From measurements of the Ca2+ dependence of PO at very negative potentials, the authors computed a variable they called RO that is equivalent to ZOR/ZCR, and by plotting logRO versus [Ca2+] obtained log〈C〉. They also obtained the intrinsic pore charge ΔqL by fitting the negative asymptote of log(PO). They recognized the importance of using log(PO), or the floor of the G-V curve, to adequately constrain parameters at negative potentials far from the midpoint of activation, although their method for estimating PO (using an all-points histogram of single-channel currents to measure NPO and dividing by NGmax/gL) did not lend itself to accurate measurements at the ceiling of the G-V curve (a requirement for defining the positive V asymptote of the conductance Hill plot). They cited this as a difficulty in estimating the value of D, defined in linkage terms as the height difference in the two asymptotes of WH[g] versus V. Instead, an alternative method for estimating D using the product of ΔqJ and the midpoint voltages of 〈qC〉 and 〈qO〉 was proposed specifically for the HA model (Scheme 2). The proof for the more general Scheme 3 goes as follows: recognizing that 〈qO〉 − 〈qC〉 = dWH[L]/dV, we integrate both sides with respect to V to obtain ∫〈qOdV − ∫〈qCdV = VΔWH[g]. Separating 〈qO〉 into its baseline (ΔqL) and variable (〈qO〉 − ΔqL) components, we subtract VΔμL = ΔqL(V(+)V(−)) from both sides of the equation and integrate the remaining expression on the left side by parts, obtaining −ΔqJ(VM(O)VM(C)) = VΔ(WH[g] + ηL). VM(C) and VM(O) are the first moments of normalized capacities derived from 〈qC〉 and 〈qO〉 − ΔqL, respectively. Recalling that VΔ(WH[g] + ηL) = −4kTln(D), and defining LΔVM = VM(O)VM(C), we obtain the desired result: ΔqJ LΔVM = 4ΔWD (see Fig. 15). Assuming one can measure both 〈qC〉 and 〈qO〉—for example, by integrating the fast component of gating charge after prepulsing to either V(−) or V(+)—then ΔWD could be calculated without resorting to logarithmic (Hill) analysis of the G-V curve. The value of ΔqJ, which is necessary to complete the calculation, can theoretically be obtained by subtracting ΔqL from ΔqT, both of which can be independently measured. The authors acquired 〈qO〉 indirectly through the relation 〈qO〉 = 〈qa〉 + 〈q〉, where 〈qa〉 = kTdlnG/dV is sensitive to the floor of the G-V curve. They acknowledged that error from this measurement could account for the 13% difference in log(D) compared with the fitted value from log(PO)-V curves.

The allosteric factor E was estimated in similar fashion by measuring shifts in 〈qC〉 and 〈qO〉 in response to saturating changes in [Ca2+]i. Because the influence of pore activation has been removed from both 〈qC〉 and 〈qO〉, the two curves can be described by a Boltzmann function in the HA model (and very nearly so in Scheme 4 because of the very large negative value of ΔWG coupling the two J particles; see Fig. 15). The success of the HA model in fitting ionic and gating current data, particularly at voltages where PO was very low (∼10−7), is reassuring for the eventual success of linkage methods described here. It will be of interest to see if, despite its potential limitations, the slow voltage ramp can overcome the inability of the NPO histogram to generate positive voltage (“ceiling”) asymptotes in the conductance Hill plot.

### Deconstructing other thermodynamic methods of analysis

Methods for dealing with equilibrium plots are prevalent in the ion channel literature. The tools used most often are the Hill and Boltzmann equations, which can provide very good fits to experimental data. Confusion arises from different uses of the “Hill” label. The traditional Hill plot is defined as ln(M/(MmaxM)), where M is a marker of activation for the principal protomer (in titration experiments, for which Hill analysis was developed, M is the binding curve). The source of differences lies in the models to which Hill analysis is applied. In the present work, allosteric interaction between the principal component (usually the pore, marked by the conductance G) and other particles generates a sigmoidal curve whose linear asymptotes reflect the strength of the interaction. A nearly identical treatment is the χ-analysis of Chowdhury and Chanda (2010), which was proposed as a means of interpreting mutational effects on local interactions between the principal particle and its neighbors.

On the other hand, Hill (1910) derived an empirical formula to describe ligand titration curves based on the partition function Z = 1 + Kn, which, interpreted in molecular terms, implies concerted movement among n protomers. This is an idealized construct that occasionally approaches reality in the setting of strong homologous interactions (for example, open–closed transitions in the tetrameric ion channel pore, as explained earlier). In the language of this paper, the empirical formula describes a single particle X whose equilibrium constant is X = Kn and whose activation curve is therefore Kn/(1 + Kn), known as the “Hill equation.” The difficulty arises when the empirical formula is applied to proteins whose behaviors do not fit the definition of a “single-particle” system. Nevertheless, this is commonly done in fits of binding curves from allosteric proteins, in which n is the fitted parameter, named nH or the “Hill coefficient,” and is used as an estimate for the minimum number of binding sites with respect to the principal ligand. The Hill equation, whose “Hill plot” is a straight line, effectively ignores the asymptotic behavior of Hill-transformed data and focuses on the maximum slope of the rapidly rising sigmoid-shaped region. In voltage-sensitive systems, such as the tetrameric pore described earlier, fits of Q-V curves to a Boltzmann function are actually fits to the Hill equation in disguise, in which the fitted Boltzmann charge ΔqB = nHΔq, where Δq is the subunit gating charge displacement (Yifrach, 2004). This is fairly imprecise because the Q-V curve is not meant to be analyzed through Hill transformation, as it generally receives contributions from more than one species of charged particles. For any distribution of charged particles more complex than a group of identical copies, The Q-V curve is more properly analyzed through global capacitance analysis using the midpoint parameter VM and the expression WC[q] = −ΔqTVM, as described here and in Chowdhury and Chanda (2012).

Another conventional use of the “Boltzmann technique” is to detect curve shifts on the V axis in response to a perturbation. The “half-activation” voltage V1/2, as determined by the midpoint of a Boltzmann curve fitted to the data, is sometimes used as a measure of the effort required to activate a channel. In analyzing Q-V curves, V1/2 is often used as a stand-in for the more correct VM. In a somewhat confusing twist of nomenclature, the V1/2 (the voltage for which the normalized Q-V is 0.5) is equal to the median value of the capacitance distribution function fq (see accompanying article by Chowdhury and Chanda in this issue), whereas VM, which according to the convention established by Wyman should be called the “median” voltage of activation, is, as explained previously, its first moment, or mean value. In any case, V1/2 has no special thermodynamic significance (in the sense of describing an energy or displacement), but it may coincide with VM for symmetrical Q-V plots.

Similarly, the product ΔqBV1/2 derived from the G-V curve does not generally represent the work of opening the channel. Nevertheless, an exception of sorts can be made for Kv channels that possess a single late-activating open state, such as the well-characterized Shaker channel (Hoshi et al., 1994; Schoppa and Sigworth, 1998). The reason is that, in these channels, the gating charge can be written as 〈q〉 = ΔqT − 〈qa〉, where 〈qa〉 = kTdlnG/dV is the mean activation charge displacement. Therefore, gating current measurements are not required to measure WC[q], because even ΔqT can be obtained from the limiting value of activation charge displacement: ΔqT = 〈qaV(−). However, the logarithmic dependence of 〈qa〉 on G implies that the extreme edges or “tails” of the G-V are important in evaluating WC[q], which for a late-activating Kv channel is related to G through −kTV(d2lnG/dV2)dV. This contrasts with the traditional Boltzmann method of analysis, where ΔqBV1/2 obtained by Boltzmann fitting is numerically equal to −kTln(B/(1B))V=0. Boltzmann fits are not typically weighted toward the tail regions but, as implied earlier, are judged by their ability to determine the position and steepness of the more accessible midpoint part of the curve. Despite its significant limitations, the Boltzmann technique been shown to be useful in studying double mutant cycles involving pore residues of late-activating Kv channels (Zandany et al., 2008), for which ΔqBΔV1/2 was found to be proportional to the true perturbation energy of a pore mutation over a large range in L (Yifrach and MacKinnon, 2002).

Finally, just as it is incorrect to apply local Hill analysis to Q-V curves, it is equally wrong, though tempting, to apply global analysis to G-V curves by proposing a new linkage function WC[g] = ΔqLVdPO = ΔqLVM[g], where VM[g] = ∫VfgdV and fg = dPO/dV. The variable VM[g] could be referred to as the “mean” activation voltage of conduction, and fg would be the “conductance capacitance,” even though this is a nonsensical term. It is nevertheless reasonable, in the context of the BK channel, to enquire whether ΔqLμΔVM[g], as a kind of “global analysis” of the local K–L interaction, is equal to −4ΔWC. This attempt at deriving the strength of K–L linkage (or any perturbation affecting pore opening, say from a mutation) by measuring shifts in VM[g] is flawed, because, unlike Q = NkT(∂lnZ/∂lnV), the conductance G = (∂lnZ/∂lnL)Gmax cannot be expressed as the voltage derivative of Z (unless the pore is the only source of voltage dependence in the channel, as in Scheme 1). Attempts at deriving a G-V–based version of Eq. 32 by integrating the putative work function WC[g] by parts to obtain ΔqLPOdV − ΔqLV(+), and then using the chain rule on the first term, as in:

$ΔqL∫−∞∞(∂lnZ∂lnL)dV=ΔqL∫−∞∞(∂lnZ∂V)(∂V∂lnL)dV+ΔqL∫−∞∞(∂lnZ∂μ)(∂μ∂lnL)dV=kTln(ZV(+)ZV(−))+kTΔqLΔnL∫−∞∞(∂lnZ∂μ)dV,$

are problematic, particularly as the second term diverges for ΔnL = 0. Recognizing that PO = −kT(∂lnZ/∂ηL) points to the underlying problem, which is that ηL cannot be varied independently of ηJ because both particle potentials are functions of the control parameter V. We contrast this to the earlier correct method (using Q-V curves) of applying global capacitance analysis to the local J–K interaction to obtain −ΔWE (Eq. 35), in which the only other source of charge movement besides the J particle was eliminated by locking the pore in the closed or open position.

To summarize this section, we have excluded certain candidates from a list of valid work functions. These include: (a) the product ΔqBV1/2 derived from Boltzmann fits of both Q-V and G-V curves; (b) local (Hill) analysis of the Q-V curve; and (c) global (capacitance) analysis of the G-V curve. Exceptions for certain constrained systems are as noted above.

### Other applications of linkage analysis: Temperature gating in TRP

Linkage analysis can in principle be applied to any ion channel with an “accessible”-state space, whose regulatory particles are energetically but not obligatorily coupled. Gating schemes like those applied to BK in this paper can be generalized for use with other tetrameric channels simply by adjusting the system equation (Eq. 3) and particle potentials (Eq. 6) to reflect the mix of external forces in play. Of special interest are temperature-regulated channels, of which there are several varieties in the TRP family (Latorre et al., 2009). Although all proteins are to some degree temperature dependent as a consequence of being immersed in a heat bath, the ability to activate a channel within a very narrow range in temperature might require a special sensor. There is considerable uncertainty about where such a T-sensor might be located (Latorre et al., 2009; Yao et al., 2011; Cui et al., 2012), and whether there is more than one site of temperature sensitivity (Clapham and Miller, 2011). Assuming that a particular domain K can be identified as a T-sensor, its particle potential will have the minimum form: ηK = ΔHK − ΔSKT, where ΔH and ΔS are the activation enthalpy and entropy, respectively, of a mole of K particles, and its activation, or “melting,” curve will be given by K/(1 + K), where the equilibrium constant K = exp(−ηK/RT). The midpoint of the curve, or “melting temperature,” is given by To = ΔHKSK. The steepness of the curve is determined by ΔSK. Rearranging K to read

$K=RΩ(T−T0T),$

we can estimate the proportional increase in the density of states, RΩ = exp(ΔSK/R), required to “melt” the T-sensor (arbitrarily defined as a change in K/(1 + K) from 0.1 to 0.9) in response to a 10°C change in temperature, a sensitivity that is achievable by both “cold”-activated (Brauchi et al., 2004) and “hot”-activated (Yao et al., 2010) TRP channels. This turns out to be ∼1055, which is equivalent to an entropy change of ΔSK = 251 cal/mol/K. Such a large increase in the total configurational space (protein plus lipids plus solution) is most compatible with a denaturing process. It has been estimated that protein folding requires a configurational entropy change of about −2.9 cal/mol/K per residue (Pickett and Sternberg, 1993). This suggests that an 87-residue peptide, with a folding enthalpy ΔHK = ΔSKTo ≈ −74 kcal/mol, could account for the steepness of the melting curve. When divided between four subunits, this would be ∼22 residues per subunit, a reasonable number for a protein domain.

Assuming that the T-sensor is allosterically linked to pore opening—perhaps through a model such as Scheme 1, or if including voltage sensors, Scheme 2—then linkage techniques can be used to determine the strength of the coupling energy ΔWC. The enthalpic component of ΔWC can be measured from the rise of the sigmoidal component of the conductance Hill energy WH[g] plotted versus T. If there is also an entropic component ΔSC, it would show as a difference in the slopes of the positive and negative asymptotes equal to ΔSC.

A distinguishing feature that sets a putative T-sensor apart from the “typical” voltage sensor is the substantial increase in heat capacitance of the denatured peptide compared with the native form. Recall from the earlier discussion that the state-dependent “fast” gating charge capacitance in Shaker is fairly insignificant compared with peak gating capacitance (∼0.4%), leading to only very subtle differences in the limiting slopes of the Q-V curve (Sigg et al., 2003). In contrast, the increase in heat capacitance ΔCp experienced by a protein upon denaturing is typically about five times the value of ΔSKo evaluated at To (Privalov and Khechinashvili, 1974), which would generate a measurable nonlinearity in the particle potential ηK of a T-sensor. The value of ΔCp appears to be fairly constant (Baldwin, 1986), leading to ηK = (ΔCp − ΔSKo)(TTo) − ΔCpTln(T/To), which is derived by integrating the thermodynamic identities Cp = (∂H/∂T)P = T(∂S/∂T)P around To = ΔHKoSKo, evaluated for both native and denatured states (Clapham and Miller, 2011).

As a consequence of the change in heat capacitance with melting, dramatic changes are seen in WH[g] related to the possibility that ηK crosses zero twice, theoretically generating two melting temperatures separated approximately by 2ΔT = 2ΔHKoCp, where ΔT is the difference between the zero energy and zero entropy temperatures of the particle potential (Schellman, 1987; see derivation in the supplemental text). At the larger temperature To, the sensor would be hot-activated (negative slope for ηK), whereas at the lower temperature, it would be cold-activated, suggesting that both phenotypes could theoretically exist in the same channel (Clapham and Miller, 2011). This is one of three explanations for phenotype reversal in domain-swapping experiments (Brauchi et al., 2006) examined in the supplemental text as an illustration of how linkage analysis might discriminate between mechanisms for temperature gating. Because of the conductance Hill plot’s logarithmic dependence on temperature, it is well suited to detect subtle differences in the shapes of a G-T curve (particularly in the asymptotic regions), as compared with the traditional Boltzmann curve–fitting technique.

We neglect to discuss here the small but real voltage dependence of TRP gating, which is dealt with elsewhere (Matta and Ahern, 2007; Latorre et al., 2009). Needless to say, it is straightforward to insert a −ΔqV term in the particle potential of any regulatory particle, whether as part of a separate voltage sensor, such as the J particle in Schemes 2–5, or to infuse voltage sensitivity into the L and K particles.

### Comparison of thermodynamic- and kinetic-based models and methods

Beginning with the well-known Hodgkin and Huxley equations describing action potentials in the squid giant axon (Hodgkin and Huxley, 1952), the interpretation of electrophysiological data has traditionally been done with the help of kinetic models. There are several reasons for this. The telegraph-like signal of single-channel ionic recordings and the multi-exponential transients of gating currents invite a kinetic description. Some ion channels enter into inactivated states when subjected to extended membrane depolarization, and kinetic analysis can be used to differentiate between the different modes of gating. Finally, many channels appear to have a fairly small number of kinetically distinguishable states, making it easier to construct and choose among competing kinetic models.

On the other hand, there are advantages to considering thermodynamic models and the linkage techniques used to analyze them. The measured quantities in linkage analysis are energy and displacement, which are used to characterize gating networks consisting of regulatory particles. Kinetic analysis, on the other hand, measures rate constants and displacements for a particular scheme of connected microscopic states. The partition function, which defines the thermodynamic properties of a channel, can often be expressed as a simple polynomial, even for complex allosteric models such as Schemes 4 and 5 that possess >103 states, whereas the equivalent kinetic models in these schemes are cumbersome to create and analyze (consider the kinetic version of Scheme 4 described in Materials and methods).

In thermodynamic modeling, coarse-graining the partition function to match the state space with the available data is straightforward; one simply sums over states that cannot be directly observed. The resulting CPFs play a central role in evaluating linkage functions, as discussed in this paper. In kinetic models, coarse-graining is effectively performed by resolving time constants of decay, which, because they are usually composed of multiple elementary processes, may be difficult to ascribe to a particular structure or mechanism of action. To do so requires at least one of the following: (a) the observation of a kinetically distinct process through sheer luck; (b) the analysis of single-channel data or measuring fluctuations; or (c) recording under extreme conditions to achieve a “limiting rate.” The last is technically difficult to do in a step-pulse experiment but constitutes precisely the asymptotic region where linkage analysis operates (and where, coincidentally, there is the least “lag time” in ramp experiments). With some exceptions (Sigg et al., 2003; Chakrapani and Auerbach, 2005), kinetic analysis operates in the configurational “middle ground” where time constants are slowest. In this respect, the two approaches complement each other.

The number of independent variables in a thermodynamic model is much smaller than in the corresponding kinetic model. A simple voltage-dependent transition requires only two thermodynamic parameters (ΔG and Δq) but four kinetic parameters (α, β, Δqα, and Δqβ). In allosteric systems, interaction energies are simply added to existing particle potentials, but in kinetic models, one must decide how to apportion the interaction energy between the transition rate constants α and β. This is typically done with respect to the position of the transition barrier in accordance with the theory of linear free energy relationships (Fersht et al., 1987; Auerbach, 2005). Using linear free energy relationship methods, it has been demonstrated that the transition state of pore activation in Shaker closely resembles the open state (Azaria et al., 2010). This suggests that the distribution of open times obtained from single-channel analysis could be wholly insensitive to allosteric perturbations of pore opening, which would instead be reflected in the more complex distribution of closed times. In linkage analysis, the position of the transition state is irrelevant, provided adequate time is allowed for equilibration. To be fair, thermodynamic methods have, as discussed earlier, failed in the case of Shaker to measure the interaction energy between pore and voltage sensor, as inferred, for example, from the inability to measure a reduction in 〈qa〉 at very negative membrane potentials (Islas and Sigworth, 1999). This outcome is not surprising from a kinetic standpoint, as single-channel analysis of Shaker reveals only one open state and many closed states (Hoshi et al., 1994), implying, from an allosteric standpoint, that there are many open states in Shaker that are energetically inaccessible.

The preceding example stresses the importance of kinetic analysis in developing models of gating, particularly in channels like Shaker with a relatively constrained-state space. However, in channels with a larger-state space, linkage analysis can characterize gating networks in a manner that speaks to the energetic interactions between activated states (particles) in regulatory domains. This is substantively different from the traditional goal of kinetic modeling in ion channels, which is to construct a rate constant diagram of accessible microscopic states. As new markers of activation are developed and refined, we may be poised to see an expansion of the application of rigorous thermodynamic methods to study energy relationships between experimentally identified sites of environmental sensitivity and the pore.

Linkage diagrams are graphical descriptions of allosteric models developed for this paper and are to my knowledge not used elsewhere. They are meant to be read from left to right. The arrow symbol (→) represents a change in the equilibrium constant of one particle mediated by the activation of another particle, whose equilibrium constant is located above the arrow. Linkage diagrams have several useful features. First, the diagrams are precise (although not unique) representations of allosteric models, leaving no ambiguity as to the arrangement of particle interactions except when bold-type variables are used to represent complex domains. Second, they are intimately connected to the partition function Z and in fact recapitulate the parsing procedure used to write down Z in polynomial form, as determined by the set of rules outlined in Table 3. Third, they are easily generated using equation editors and take up little vertical space on the page.

In linkage notation, CPFs are represented by quantities within parentheses. The simplest CPF is that of a single particle K: (K)⇒1 + K. Multiplying the equilibrium constant K by the allosteric factor C yields (KC)⇒1 + KC. Mediating this change in K might be the particle L, leading to $(K)→L(KC)$⇒(1 + K) + L(1 + KC). If C = 1, implying independence between L and K, then the CPF becomes (1 + K)(1 + L), which can also be written in linkage notation as (K)(L). Using these rules and the linkage diagrams found in the text, the partition function for Scheme 1 is seen to be Z = (1 + K)4 + L(1 + KC)4, and for Scheme 2, it is Z = [(1 + K) + J(1 + KE)]4 + L[(1 + KC) + JD(1 + KCE)]4. Linkage diagrams are not unique. For example, $(K)→L (KC)$ is equivalent to $(L)→K (LC)$. For an allosteric model with m particle species, there are m! ways of drawing its linkage diagram and, by extension, m! ways of writing down its partition function. Complex domains consisting of more than one particle can be specified through the use of bold type. For example, a gating ring domain with multiple Ca2+ sensors K1, K2, K3… can be represented by K[F], where F = {F1, F2, F3…} are internal allosteric factors. In cases where the internal allosteric network is not known, or the domain is otherwise complex, it is not possible to write down the corresponding CPF explicitly. However, the two extreme terms of the CPF may still be evaluated. For K[F], they are (K[F])(−) = 1 and (K[F])(+) = (K1K2K3F1F2F3…). Because linkage analysis is chiefly concerned with these extreme values, the linkage diagram of complex systems remains a useful tool even if the partition function cannot be fully defined.

I am grateful to my friend and longtime collaborator Riccardo Olcese, and members of his laboratory, for turning me on to the subject of BK channels, and for helpful comments on the manuscript. Heartfelt thanks also to Baron Chanda and Sandipan Chowdhury, who wrote the companion piece to this paper, for reviewing the revised manuscript, and to two anonymous reviewers and the editorial staff of the JGP.

Christopher Miller served as editor.

Aggarwal
S.K.
,
MacKinnon
R.
.
1996
.
Contribution of the S4 segment to gating charge in the Shaker K+ channel
.
Neuron.
16
:
1169
1177
.
Alberty
R.A.
2002
.
Use of legendre transforms in chemical thermodynamics: International Union of Pure and Applied Chemistry, Physical Chemistry Division, Commission on Thermodynamics
.
J. Chem. Thermodyn.
34
:
1787
1823
.
Auerbach
A.
2005
.
Gating of acetylcholine receptor channels: Brownian motion across a broad transition state
.
102
:
1408
1412
.
Azaria
R.
,
Irit
O.
,
Ben-Abu
Y.
,
Yifrach
O.
.
2010
.
Probing the transition state of the allosteric pathway of the Shaker Kv channel pore by linear free-energy relations
.
J. Mol. Biol.
403
:
167
173
.
Baldwin
R.L.
1986
.
Temperature dependence of the hydrophobic interaction in protein folding
.
83
:
8069
8072
.
Brauchi
S.
,
Orio
P.
,
Latorre
R.
.
2004
.
Clues to understanding cold sensation: thermodynamics and electrophysiological analysis of the cold receptor TRPM8
.
101
:
15494
15499
.
Brauchi
S.
,
Orta
G.
,
Salazar
M.
,
Rosenmann
E.
,
Latorre
R.
.
2006
.
A hot-sensing cold receptor: C-terminal domain determines thermosensation in transient receptor potential channels
.
J. Neurosci.
26
:
4835
4840
.
Chakrapani
S.
,
Auerbach
A.
.
2005
.
A speed limit for conformational change of an allosteric membrane protein
.
102
:
87
92
.
Chapman
M.L.
,
VanDongen
H.M.
,
VanDongen
A.M.
.
1997
.
Activation-dependent subconductance levels in the drk1 K channel suggest a subunit basis for ion permeation and gating
.
Biophys. J.
72
:
708
719
.
Chowdhury
S.
,
Chanda
B.
.
2010
.
Deconstructing thermodynamic parameters of a coupled system from site-specific observables
.
107
:
18856
18861
.
Chowdhury
S.
,
Chanda
B.
.
2012
.
Estimating the voltage-dependent free energy change of ion channels using the median voltage for activation
.
J. Gen. Physiol.
139
:
3
17
.
Clapham
D.E.
,
Miller
C.
.
2011
.
A thermodynamic framework for understanding temperature sensing by transient receptor potential (TRP) channels
.
108
:
19492
19497
.
Conti
F.
,
1986
.
The relationship between electrophysiological data and thermodynamics of ion channel conformations
.
In
Ion Channels in Neural Membranes.
Bolis
L.
,
Keynes
R.D.
,
editors
.
Alan R. Liss, Inc
,
New York
.
25
41
.
Cox
D.H.
,
Cui
J.
,
Aldrich
R.W.
.
1997
.
Allosteric gating of a large conductance Ca-activated K+ channel
.
J. Gen. Physiol.
110
:
257
281
.
Crouzy
S.C.
,
Sigworth
F.J.
.
1993
.
Fluctuations in ion channel gating currents. Analysis of nonstationary shot noise
.
Biophys. J.
64
:
68
76
.
Cui
Y.
,
Yang
F.
,
Cao
V.
,
Yarov-Yarovoy
V.
,
Wang
K.
,
Zheng
J.
.
2012
.
Selective disruption of high sensitivity heat activation but not capsaicin activation of TRPV1 channels by pore turret mutations
.
J. Gen. Physiol.
139
:
273
283
.
De Wet
H.
,
Allen
M.
,
Holmes
C.
,
Stobbart
M.
,
Lippiat
J.D.
,
Callaghan
R.
.
2006
.
Modulation of the BK channel by estrogens: examination at single channel level
.
Mol. Membr. Biol.
23
:
420
429
.
Di Cera
E.
1990
.
Thermodynamics of local linkage effects. Contracted partition functions and the analysis of site-specific energetics
.
Biophys. Chem.
37
:
147
164
.
Di Cera
E.
,
Chen
Z.Q.
.
1993
.
The binding capacity is a probability density function
.
Biophys. J.
65
:
164
170
.
Di Cera
E.
,
Gill
S.J.
,
Wyman
J.
.
1988
.
Binding capacity: cooperativity and buffering in biopolymers
.
85
:
449
452
.
Eaton
W.A.
,
Henry
E.R.
,
Hofrichter
J.
,
Bettati
S.
,
Viappiani
C.
,
Mozzarelli
A.
.
2007
.
Evolution of allosteric models for hemoglobin
.
IUBMB Life.
59
:
586
599
.
Ferguson
W.B.
,
McManus
O.B.
,
Magleby
K.L.
.
1993
.
Opening and closing transitions for BK channels often occur in two steps via sojourns through a brief lifetime subconductance state
.
Biophys. J.
65
:
702
714
.
Fernández
J.M.
,
Bezanilla
F.
,
Taylor
R.E.
.
1982
.
Distribution and kinetics of membrane dielectric polarization. II. Frequency domain studies of gating currents
.
J. Gen. Physiol.
79
:
41
67
.
Fersht
A.R.
,
Leatherbarrow
R.J.
,
Wells
T.N.
.
1987
.
Structure-activity relationships in engineered proteins: analysis of use of binding energy by linear free energy relationships
.
Biochemistry.
26
:
6030
6038
.
Gagnon
D.G.
,
Bezanilla
F.
.
2009
.
A single charged voltage sensor is capable of gating the Shaker K+ channel
.
J. Gen. Physiol.
133
:
467
483
.
Gamal El-Din
T.M.
,
Heldstab
H.
,
Lehmann
C.
,
Greef
N.G.
.
2010
.
Double gaps along Shaker S4 demonstrate omega currents at three different closed states
.
Channels (Austin).
4
:
93
100
.
Henrion
U.
,
Renhorn
J.
,
Börjesson
S.I.
,
Nelson
E.M.
,
Schwaiger
C.S.
,
Bjelkmar
P.
,
Wallner
B.
,
Lindahl
E.
,
Elinder
F.
.
2012
.
Tracking a complete voltage-sensor cycle with metal-ion bridges
.
109
:
8552
8557
.
Hill
A.V.
1910
.
The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves
.
J. Physiol.
40
:
iv
vii
.
Hille
B.
2001
.
Ionic Channels in Excitable Membranes. Third edition
.
Sinauer Associates, Inc.
,
Sunderland, MA
.
814 pp
.
Hodgkin
A.L.
,
Huxley
A.F.
.
1952
.
A quantitative description of membrane current and its application to conduction and excitation in nerve
.
J. Physiol.
117
:
500
544
.
Horn
R.
,
Ding
S.
,
Gruber
H.J.
.
2000
.
Immobilizing the moving parts of voltage-gated ion channels
.
J. Gen. Physiol.
116
:
461
476
.
Horrigan
F.T.
,
Aldrich
R.W.
.
2002
.
Coupling between voltage sensor activation, Ca2+ binding and channel opening in large conductance (BK) potassium channels
.
J. Gen. Physiol.
120
:
267
305
.
Horrigan
F.T.
,
Ma
Z.
.
2008
.
Mg2+ enhances voltage sensor/gate coupling in BK channels
.
J. Gen. Physiol.
131
:
13
32
.
Hoshi
T.
,
Zagotta
W.N.
,
Aldrich
R.W.
.
1994
.
Shaker potassium channel gating. I: Transitions near the open state
.
J. Gen. Physiol.
103
:
249
278
.
Islas
L.D.
,
Sigworth
F.J.
.
1999
.
Voltage sensitivity and gating charge in Shaker and Shab family potassium channels
.
J. Gen. Physiol.
114
:
723
742
.
Koshland
D.E.
Jr
,
Némethy
G.
,
Filmer
D.
.
1966
.
Comparison of experimental binding data and theoretical models in proteins containing subunits
.
Biochemistry.
5
:
365
385
.
Latorre
R.
,
Zaelzer
C.
,
Brauchi
S.
.
2009
.
Structure-functional intimacies of transient receptor potential channels
.
Q. Rev. Biophys.
42
:
201
246
.
Ledwell
J.L.
,
Aldrich
R.W.
.
1999
.
Mutations in the S4 region isolate the final voltage-dependent cooperative step in potassium channel activation
.
J. Gen. Physiol.
113
:
389
414
.
Long
S.B.
,
Campbell
E.B.
,
Mackinnon
R.
.
2005
.
Voltage sensor of Kv1.2: structural basis of electromechanical coupling
.
Science.
309
:
903
908
.
Markin
V.S.
,
Sachs
F.
.
2004
.
Thermodynamics of mechanosensitivity
.
Phys. Biol.
1
:
110
124
.
Matta
J.A.
,
Ahern
G.P.
.
2007
.
Voltage is a partial activator of rat thermosensitive TRP channels
.
J. Physiol.
585
:
469
482
.
Mistry
D.K.
,
Garland
C.J.
.
1998
.
Characteristics of single, large-conductance calcium-dependent potassium channels (BKCa) from smooth muscle cells isolated from the rabbit mesenteric artery
.
J. Membr. Biol.
164
:
125
138
.
Moczydlowski
E.
,
Alvarez
O.
,
Vergara
C.
,
Latorre
R.
.
1985
.
Effect of phospholipid surface charge on the conductance and gating of a Ca2+-activated K+ channel in planar lipid bilayers
.
J. Membr. Biol.
83
:
273
282
.
Monod
J.
,
Wyman
J.
,
Changeux
J.P.
.
1965
.
On the nature of allosteric transitions: a plausible model
.
J. Mol. Biol.
12
:
88
118
.
Nagel
G.
,
Ollig
D.
,
Fuhrmann
M.
,
Kateriya
S.
,
Musti
A.M.
,
Bamberg
E.
,
Hegemann
P.
.
2002
.
Channelrhodopsin-1: a light-gated proton channel in green algae
.
Science.
296
:
2395
2398
.
Pantazis
A.
,
Gudzenko
V.
,
Savalli
N.
,
Sigg
D.
,
Olcese
R.
.
2010
.
Operation of the voltage sensor of a human voltage- and Ca2+-activated K+ channel
.
107
:
4459
4464
.
Pickett
S.D.
,
Sternberg
M.J.
.
1993
.
Empirical scale of side-chain conformational entropy in protein folding
.
J. Mol. Biol.
231
:
825
839
.
Piskorowski
R.
,
Aldrich
R.W.
.
2002
.
Calcium activation of BK(Ca) potassium channels lacking the calcium bowl and RCK domains
.
Nature.
420
:
499
502
.
Press
W.
,
Teukolsky
S.
,
Vetterling
W.
,
Flannery
B.
.
1992
.
Numerical Recipes in C. Second edition
.
Cambridge University Press
,
New York
.
994 pp
.
Privalov
P.L.
,
Khechinashvili
N.N.
.
1974
.
A thermodynamic approach to the problem of stabilization of globular protein structure: a calorimetric study
.
J. Mol. Biol.
86
:
665
684
.
Qian
X.
,
Niu
X.
,
Magleby
K.L.
.
2006
.
Intra- and intersubunit cooperativity in activation of BK channels by Ca2+
.
J. Gen. Physiol.
128
:
389
404
.
Rothberg
B.S.
,
Magleby
K.L.
.
1999
.
Gating kinetics of single large-conductance Ca2+-activated K+ channels in high Ca2+ suggest a two-tiered allosteric gating mechanism
.
J. Gen. Physiol.
114
:
93
124
.
Rothberg
B.S.
,
Magleby
K.L.
.
2000
.
Voltage and Ca2+ activation of single large-conductance Ca2+-activated K+ channels described by a two-tiered allosteric gating mechanism
.
J. Gen. Physiol.
116
:
75
99
.
E.
,
Yifrach
O.
.
2007
.
Principles underlying energetic coupling along an allosteric communication trajectory of a voltage-activated K+ channel
.
104
:
19813
19818
.
Savalli
N.
,
Pantazis
A.
,
Yusifov
T.
,
Sigg
D.
,
Olcese
R.
.
2012
.
The contribution of RCK domains to human BK channel allosteric activation
.
J. Biol. Chem.
287
:
21741
21750
.
Schellman
J.A.
1987
.
The thermodynamic stability of proteins
.
Annu. Rev. Biophys. Biophys. Chem.
16
:
115
137
.
Schoppa
N.E.
,
Sigworth
F.J.
.
1998
.
Activation of Shaker potassium channels. III. An activation gating model for wild-type and V2 mutant channels
.
J. Gen. Physiol.
111
:
313
342
.
Schwappach
B.
2008
.
An overview of trafficking and assembly of neurotransmitter receptors and ion channels (Review)
.
Mol. Membr. Biol.
25
:
270
278
.
Seoh
S.A.
,
Sigg
D.
,
Papazian
D.M.
,
Bezanilla
F.
.
1996
.
Voltage-sensing residues in the S2 and S4 segments of the Shaker K+ channel
.
Neuron.
16
:
1159
1167
.
Shi
J.
,
Krishnamoorthy
G.
,
Yang
Y.
,
Hu
L.
,
Chaturvedi
N.
,
Harilal
D.
,
Qin
J.
,
Cui
J.
.
2002
.
Mechanism of magnesium activation of calcium-activated potassium channels
.
Nature.
418
:
876
880
.
Sigg
D.
,
Bezanilla
F.
.
1997
.
Total charge movement per channel. The relation between gating charge displacement and the voltage sensitivity of activation
.
J. Gen. Physiol.
109
:
27
39
.
Sigg
D.
,
Qian
H.
,
Bezanilla
F.
.
1999
.
Kramers’ diffusion theory applied to gating kinetics of voltage-dependent ion channels
.
Biophys. J.
76
:
782
803
.
Sigg
D.
,
Bezanilla
F.
,
Stefani
E.
.
2003
.
Fast gating in the Shaker K+ channel and the energy landscape of activation
.
100
:
7611
7615
.
Sigworth
F.J.
1980
.
The variance of sodium current fluctuations at the node of Ranvier
.
J. Physiol.
307
:
97
129
.
Starace
D.M.
,
Bezanilla
F.
.
2004
.
A proton pore in a potassium channel voltage sensor reveals a focused electric field
.
Nature.
427
:
548
553
.
Stockbridge
L.L.
,
French
A.S.
,
Man
S.F.
.
1991
.
Subconductance states in calcium-activated potassium channels from canine airway smooth muscle
.
Biochim. Biophys. Acta.
1064
:
212
218
.
Sweet
T.B.
,
Cox
D.H.
.
2008
.
Measurements of the BKCa channel’s high-affinity Ca2+ binding constants: Effects of membrane voltage
.
J. Gen. Physiol.
132
:
491
505
.
Taglialatela
M.
,
Toro
L.
,
Stefani
E.
.
1992
.
Novel voltage clamp to record small, fast currents from ion channels expressed in Xenopus oocytes
.
Biophys. J.
61
:
78
82
.
Tombola
F.
,
Pathak
M.M.
,
Isacoff
E.Y.
.
2005
.
Voltage-sensing arginines in a potassium channel permeate and occlude cation-selective pores
.
Neuron.
45
:
379
388
.
Wang
L.
,
Sigworth
F.J.
.
2009
.
Structure of the BK potassium channel in a lipid membrane from electron cryomicroscopy
.
Nature.
461
:
292
295
.
Wyman
J.
Jr
.
1964
.
Linked functions and reciprocal effects in hemoglobin: a second look
.
19
:
223
286
.
Wyman
J.
1975
.
A group of thermodynamic potentials applicable to ligand binding by a polyfunctional macromolecule
.
72
:
1464
1468
.
Xia
X.M.
,
Zeng
X.
,
Lingle
C.J.
.
2002
.
Multiple regulatory sites in large-conductance calcium-activated potassium channels
.
Nature.
418
:
880
884
.
Yang
H.
,
Shi
J.
,
Zhang
G.
,
Yang
J.
,
Delaloye
K.
,
Cui
J.
.
2008
.
Activation of Slo1 BK channels by Mg2+ coordinated between the voltage sensor and RCK1 domains
.
Nat. Struct. Mol. Biol.
15
:
1152
1159
.
Yao
J.
,
Liu
B.
,
Qin
F.
.
2010
.
Kinetic and energetic analysis of thermally activated TRPV1 channels
.
Biophys. J.
99
:
1743
1753
.
Yao
J.
,
Liu
B.
,
Qin
F.
.
2011
.
Modular thermal sensors in temperature-gated transient receptor potential (TRP) channels
.
108
:
11109
11114
.
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.
,
MacKinnon
R.
.
2002
.
Energetics of pore opening in a voltage-gated K(+) channel
.
Cell.
111
:
231
239
.
Yuan
P.
,
Leonetti
M.D.
,
Pico
A.R.
,
Hsiung
Y.
,
MacKinnon
R.
.
2010
.
Structure of the human BK channel Ca2+-activation apparatus at 3.0 A resolution
.
Science.
329
:
182
186
.
Zandany
N.
,
M.
,
Orr
I.
,
Yifrach
O.
.
2008
.
Direct analysis of cooperativity in multisubunit allosteric proteins
.
105
:
11697
11702
.
Zeng
X.H.
,
Xia
X.M.
,
Lingle
C.J.
.
2005
.
Divalent cation sensitivity of BK channel activation supports the existence of three distinct binding sites
.
J. Gen. Physiol.
125
:
273
286
.
Zheng
J.
,
Sigworth
F.J.
.
1997
.
Selectivity changes during activation of mutant Shaker potassium channels
.
J. Gen. Physiol.
110
:
101
117
.
1

100 meV = 2.30 kcal/mol = 9.65 kJ/mol.

2

The two cases (electrical and ligand-binding) are not exactly comparable. Although it is generally assumed that binding numbers Δn all have unity value, charge displacements Δq differ from particle to particle, adding an extra level of complexity to the coefficients of the partition function.

3

The notion of a fractional Δn is difficult to reconcile with the usual notion that a ligand site is either bound or not (Δn = 1). However, one can image instances in which a weakly or nonspecifically bound ligand is only partially taken out of solution, and the change in effort to maintain a constant chemical potential is therefore also partial. I am not aware of any studies that support partial binding, but in any case, the conclusions of this paper are unaffected if it can be universally proven that Δn = 1.

Abbreviations used in this paper:

• BK

large-conductance voltage- and Ca2+-dependent

•
• CPF

contracted partition function

•
• HA

Horrigan–Aldrich

•
• LPF

limiting partition function

•
• MWC

Monod–Wyman–Changeux

•
• S/N

signal to noise

•
• TRP

transient receptor potential