A linkage analysis toolkit for studying allosteric networks in ion channels

Ion channels are pore-forming proteins that regulate (gate) the flux of selected ions (K⁺, Na⁺, Ca²⁺, 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 I_g = 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 Ca²⁺-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 Ca²⁺ 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 Ca²⁺-binding sites) generates a multitude of possible channel states, the bulk of which appear to be readily accessible through experimental application of V and [Ca²⁺]_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 Ca²⁺-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 W_XY linking the activated state of a protomer Y to that of a neighboring protomer X equals the energy of the reverse interaction W_YX. 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 W_A 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 W_A = Φ_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 W_AB and not the total energy of activating both components (= W_A + W_B + W_AB). This is helpful because, in many instances, it is not possible to measure W_B, 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 F_B 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 W_AB between components A and B can be obtained by twice measuring W_A (i.e., pulling the A lever), once with component B in the resting state (W_A,B(−)) and a second time with component B in the activated state (W_A,B(+)). The difference between these two energies yields W_AB, the logic being that if A and B are allosterically linked, activating B will assist (or retard) activation A by the amount W_AB. This form of linkage analysis can be concisely summarized by the expression W_AB = ^BΔW_A = ^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 W_AB = (Φ_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 W_A, 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 W_B and ^AΔ as its legs, should yield the same interaction energy; in other words, W_AB = W_BA. 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 W_ABC (in which the activation state of a third component C influences W_AB), are not covered here but represent a straightforward extension of the theory.

The allosteric relationship between components A and B is established by plotting W_A against F_B, referred to here as a linkage plot. The typical appearance of a linkage plot is shown in Fig. 1. The work function W_A rises steeply in the region of B activation, with a displacement equal to W_AB. If the principal component A is itself not sensitive to F_B, 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 F_B (for example, if F_B = V, then m = Δq_A, 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, W_L = Φ_L(+) − Φ_L(−). It will be shown later that W_L is derived from the G-V curve through W_L = −kTln[G/(G_max − G)]. 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, −W_L will be referred to as the “conductance Hill energy” W_H[g]. Defining W_H[g] as the negative value of −W_L 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, W_J and W_K, which are the work functions for voltage and Ca²⁺ sensor activation, respectively. The function lnε used in “χ-analysis” of mutational effects on cooperativity (Chowdhury and Chanda, 2010) is precisely equal to W_H[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 W_q of moving a channel’s gating charge q = Q/N, where Q is total measured charge, and N is the number of channels. W_q is derived by integrating Vdq over the saturating range of q given by [0, Δq_T]. Changing the integrand from q to V yields W_q = Δq_TV_M, where V_M = ∫Vf_qdV is the mean value of V with respect to the capacitance distribution f_q = (dq/dV)/Δq_T. 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 W_C[q], which we define as W_C[q] = −W_q, 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, W_n = Δn_Tμ_M describes the reversible energy of ligand binding, where μ_M is the mean value of the chemical potential μ integrated over the “binding capacity” dn/dμ (Di Cera et al., 1988), and n is the ligand number.

The summed energy W_qn of interactions between voltage-sensing and Ca²⁺-sensing elements in the BK channel is equal to the expression ^μΔW_q = Δq_T ^μΔV_M, whose evaluation in practice involves measuring the shift in V_M 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 C_q = dQ/dV, from which V_M can be obtained. Wyman’s reciprocity principle assures us that ^VΔW_n = ^μΔW_q, or more succinctly, that W_qn = W_nq. In other words, measuring the mean ligand activity μ_M as a function of V could in principle be used to confirm the value of W_qn obtained from the μ dependence of V_M, provided one is able to measure Ca²⁺-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 Ca²⁺-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 W_C[q] versus μ would demonstrate asymptotic slopes equal to Δn, the number of Ca²⁺ bound per site. However, in most models of BK, there is the assumption of strict modality separation between voltage and Ca²⁺ 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 (W_C[q],μ[−] and W_C[q],μ[+]) are required to measure W_qn = ^μΔW_q.

With the two essential work functions, W_H[g] and W_C[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.

Q-V and G-V curves for the 17-particle model (Scheme 3) were computed from particle activation curves (Eqs. 39a–c) 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 d_i based on its state-dependent transition rate α_i. During constant-voltage segments of the simulation, or if the particle charge Δq_i = 0, d_i was obtained from the formula −lnr_n/α_i, where r_n is a uniform random number drawn between 0 and 1 obtained with the long-period (>2 × 10¹⁸) random number generator of L’Ecuyer with Bays-Durham shuffle (Press et al., 1992). During the ramp phase, the d_i of charged particles was obtained from the formula (2kT/mΔq_i)ln[1−(mΔq_ilnr_n)/(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(Δq_iV/2kT) and β_oi exp(−Δq_iV/2kT), respectively. The method for determining α_i and β_i are given in the next section. The particle with the shortest d_i was allowed to transition, after which the simulation time was advanced by d_i. If between t and t + d_i there was a change in pulse protocol occurring at t₁, the simulation time reverted to t₁, 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 f_c) using a Gaussian filter and discretized (sampling frequency fs = 5f_c) in such a way that the discrete integral of the response equaled the gating charge displacement Δq_i (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 g_L and driving force (V − V_eq). The value of the equilibrium potential V_eq 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 (I_noise) was equal to the Nyquist value 4kTBg, where B is the low-pass bandwidth of a Gaussian filter (B ≈ 1.064 f_c; see Crouzy and Sigworth, 1993), and the conductance g was made up of contributions from a constant “leak” conductance g_leak, assumed to be 50.0 ± 2.5 pS (≈ 20 Gohm⁻¹), and a fluctuating term N_Og_L, where N_O was the number of open channels at a given time. Each simulated current trace (gating and ionic) contained the following terms:

where C_patch 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, V_eq, g_leak, and C_patch 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 V_r = −300 mV, ending with a tail (OFF) pulse to 0 mV.

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. Ca²⁺ binding was assumed to have the same τ_max as voltage-sensor activation. The condition of maximum τ = (α + β)⁻¹ was obtained when forward (α) and backward (β) particle activation rates were both equal to (2τ_max)⁻¹. To satisfy equilibrium and kinetic requirements, particle rate constant were computed according to α_i = (2τ_max)⁻¹(Z_i(+)/Z_i(−))^1/2 and β_i = (2τ_max)⁻¹(Z_i(−)/Z_i(+))^1/2, where Z_i(−) and Z_i(+) are the limiting channel partition functions corresponding to the resting and activated states of the i^th particle, obtained during ramp simulation from a look-up table continuously updated to account for changes in V.

Simulated experiments consisted of 100 traces each of ionic and gating currents from n = 1,000 ± 50 channels, taken at [Ca²⁺]_i = 10⁻³ µM and 10³ µM. Estimates for g_leak and C_patch were obtained by analyzing short segments (2–40 ms) of current before (I₍₋₎) and after (I₍₊₎) the beginning of the test pulse (t₁). For a square pulse of height ΔV, g_leak = (I₍₊₎ − I₍₋₎)/ΔV, where the current was averaged over each segment. In ramp experiments, initial estimates of g_leak and C_patch were simultaneously obtained by fitting I to the following function: I (t < t₁) = I₍₋₎; I (t ≥ t₁) = I₍₋₎ + mC_patch + g_leak(V(t) − V_r). Fits were obtained using the solver application in Excel (Microsoft). The values of g_leak and C_patch, in addition to V_eq (in ionic current experiments), were further refined through manual adjustment and fitting of asymptotes to achieve satisfactory shapes for both the Q-V curve = ∫(Ni_g/m)dV and the conductance Hill plot W_H[g] = kTln(G/(G_max − G), where G = I/(V − V_eq). To reduce the deleterious effects of Nyquist current noise on Q-V and W_H[g] plots, post-filtering of currents was performed after initial correction for g_leak and C_patch 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 V_M was estimated from the following integral (see Eq. 29):

Integration was performed using Riemann sums, which was adequate for the sampling rate used (δV = 0.2 mV). In W_H[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 [Ca²⁺]_i = 10⁻³ and 10³ µM. The pore gating charge Δq_L was obtained from the slope of the lower (negative) V asymptote in the high Ca²⁺ 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) = iI − I²/N (Sigworth, 1980) for both test and tail currents, with N averaged over the two estimates. The total gating charge Δq_T was obtained from the integrated positive gating current ΔQ = ∫Ni_gdV 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.

The partition function and linkage relationships for Scheme 3 (see Eqs. 37a–d) 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.

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:

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 N_b and N_f. 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}] = U − ST + PV − QV − Σ_bμ_bN_b (notation adapted from Alberty, 2002). Integrating Eq. 1 at constant values of the intensive variables, we obtain G[V,{μ_b}] = Σ_fμ_f N_f, 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}] = Nμ_BK. Following Wyman (1975), we equate the system energy Φ with μ_BK. Writing the Gibbs–Duhem equation for Eq. 1 as −SdT + VdP − QdV − Σ_bN_bdμ_b − NdΦ = 0 and solving for dΦ, we obtain:

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 Ca²⁺ 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:

In Eq. 3, the Ca²⁺ chemical potential μ = μ^ο + kTlna_Ca is a function of internal calcium activity a_Ca (in practice, assumed to be [Ca²⁺]_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,

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

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 (n_J, n_K) or 1 (n_L). 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 = J₁J₂…).

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 (η_X ≡ ^XΔΦ_X), and evaluating η_X at constant T, P, V (voltage), and μ, we obtain:

where Δq_X and Δn_X are charge and Ca²⁺ displacements of activation, and ΔG_X = Δ(U/N)_X − Δs_XT + Δv_XP 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 = ΔG_X − Δq_LV, 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:

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:

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:

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

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

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

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 + J₁ + J₂ +…; or it is composed of multiple interacting two-state particles, in which case the partition function is (for two particles) Z = (1 + J₁) + J₂(1 + J₁G), with G as a coupling factor. In both cases, the particles are connected to the equilibrium constants J₁, J₂, 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 J₂, 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 〈j₁〉 = ∂lnZ/∂lnJ₁, 〈j₂〉 = ∂lnZ/∂lnJ₂, 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 + P⁴ from which the classical Hill equation (Hill, 1910) can be derived as H = (1/4)∂lnZ/∂lnP = P⁴/(1 + P⁴). 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 = P⁴.

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.

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:

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 Z_L(−) and Z_L(+). 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 Z_C and Z_O, 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 Ca²⁺-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):

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 Z_C = (1 + K), and Z_O = (1 + KC). The ratio of the two CPFs yields the expectation value of C derived from the distribution of K particle states:

The Ca²⁺ dependence in Eq. 13 arises from the equilibrium constant K = exp(−η_K/kT) = exp[(−ΔG_K + Δn_Kμ)/kT]. Recalling that [Ca²⁺]_i = exp(μ/kT), this can be rearranged to read K = ([Ca²⁺]_i/K_d)^ΔnK, where K_d = exp(ΔG_K/Δn_KkT) is the intrinsic disassociation constant of the K-binding site. If a single ion is bound during activation, then Δn_K = 1, and we obtain the usual expression K = [Ca²⁺]_i/K_d. Because K grows monotonically with increasing μ, it is correct to substitute 〈C〉_μ for 〈C〉_K 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 −ΔW_C.

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 Hill energy is defined as the work function of a single particle. It was stated in the Introduction that the pore-specific Hill energy W_L can be derived from the G-V curve. This will now be proved. Starting with the definition W_L = Φ_L(+) − Φ_L(−), we invoke the bridge equation (Eq. 4) to obtain W_L = −kTlnZ_L(+)/Z_L(−). As discussed in the preceding section, the partition function can be written as the sum of two LPFs: Z = Z_L(−) + Z_L(+). Recognizing that each term in Z is a state probability after normalizing with Z, we write the probability P_O of channel opening as G/G_max = Z_L(+)/Z. Combining the above equations, and recalling that the conductance Hill energy was defined as W_H[g] = −W_L, we obtain, as advertised, the following Hill equation:

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

From the discussion in the previous section, we know that W_H[g] + η_L varies in sigmoid fashion from 0 to −4ΔW_C 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:

It should be noted that, because the pore in Scheme 1 is not intrinsically Ca²⁺ dependent, its particle potential η_L = ΔG_L − Δq_LV is independent of μ; i.e., ^μΔη_L = 0. Under these circumstances, η_L can be left out of Eq. 16.

Plotting W_H[g] versus V instead of μ yields a straight line given by W_H[g] = −η_L. This is the Hill plot of the pore particle in the absence of interactions with other voltage sensors. Under conditions of “zero” Ca²⁺, symbolized by μ₍₋₎, the line has slope m = Δq_L and V-intercept V_L,μ(−) = ΔG_L/Δq_L, consistent with η_L = −Δq_L(V − V_L,μ(−)). If we increase Ca²⁺ to saturating levels, then m is unchanged but the V-intercept V_L,μ(+) = (ΔG_L + 4ΔW_C)/Δq_L is left-shifted in the setting of positive allosteric interaction (ΔW_C < 0). This single-particle behavior is the source of sloping asymptotes in Hill plots.

The conductance Hill plot (Eq. 14) was derived from the equilibrium curve (G-V) of a single particle L. Suppose we could measure Ca²⁺ binding in Scheme 1 through an assay in which the signal B increments by b each time a Ca²⁺ 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:

where B = Nb〈k〉 (maximum value: B_max = 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:

(SCHEME 1, alternative form)

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

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

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

It is possible to regain from Eq. 20 the simple binary ratio of CPFs that made it possible to extract the value of −4ΔW_C 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:

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 −ΔW_C (however, notice the absence of the factor 4 compared with Eq. 16). Because the particle potential η_K = ΔG_K − Δn_Kμ is a linear function of μ, the plot of W_H[b] versus μ demonstrates sloping asymptotes (m = Δn_K), as illustrated in Fig. 1 (substituting W_H[b] for W_A, and μ for F_B). Assuming Δn_K = 1, the expressions for the μ intercepts of the two asymptotic lines are μ_K,μ(−) = ΔG_K and μ_K,μ(+) = ΔG_K + ΔW_C, from which the disassociation constant K_d is obtained through K_d = exp(μ_K/kT). We see that, stemming from linkage with the pore, there are two values for K_d (Xia et al., 2002), one for the un-liganded (free) state and the other for the bound state (K_d value left-shifted for the bound state if ΔW_C < 0).

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 Ca²⁺-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 B_k. As a result, the original Z_(k) are changed to B₀Z₍₀₎, B₀B₁Z₍₁₎, B₀B₁B₂Z₍₂₎, and so forth. The zeroth factor B₀ 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 kTlnB₁ and kTlnB₄ to Eqs. 21a and 21b, respectively. Therefore, in evaluating ^μΔW_H[b] for Scheme 1, we obtain—in addition to the K–L interaction energy −ΔW_C and the divergent term −^μΔη_K = Δn_K(μ₍₊₎ − μ₍₋₎)—a third, self-interaction term kTln(B₄/B₁), representing the difference in the work needed to activate the first and last K particles.

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

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

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 Z_C and Z_O, derived from the parenthesized contents in the linkage diagram, are now functions of J and K (Table 4). A single polynomial function f₁ describes both CPFs: Z_C = f₁(J, K, E) and Z_O = f₁(JD, KC, E). Depending on whether J or K is varied, the CPF ratio Z_O/Z_C can be interpreted as an expectation value for either D or C. In analogy to Eq. 13, we have:

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:

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 ΔW_D and ΔW_E are nonzero. Otherwise, J could be viewed as merely an accessory component of either the pore or the Ca²⁺ 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 Z_C and Z_O (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.

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

The sub-schemes for the possible configurations of voltage-sensing L and J particles are given by the four (Ca²⁺-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 Z_C and Z_O, share a single polynomial function, in this case: f₂(λ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 Z_OR/Z_OR ≈ Z_OA/Z_OA = 〈C〉_μ. Alternatively, Z_C and Z_O 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 g₂(λJ) = 1 + λJ, where λ = 1, D, E, or DE (see Table 6). Rephrasing Eq. 22c, we obtain Z_OF/Z_OF ≈ Z_OB/Z_OB = 〈D〉_V. Either parsing route, if connected to a suitable work function, also enables one to measure ΔW_E, 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 ΔW_E, 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:

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 Ca²⁺-binding elements, and is described next.

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

Integrating Eq. 26 by parts, we obtain:

Although both terms in Eq. 27 diverge, W_C[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 A₁ + A₃ in Fig. 3 A. The product Δq_TV₍₊₎ is equal to A₂ + A₃. The difference between these two areas, −(A₂ − A₁), 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 W_C[q] make use of the midpoint voltage V_M, which was defined in introductory remarks as W_q/Δq_T = −W_C[q]/Δq_T. From Eq. 26, we can write the following integral expression for V_M:

Chowdhury and Chanda (2012) recently introduced V_M 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 V_M 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 A₆ + A₇ − A₄, which represents the value of W_q/Δq_T and is therefore equal to V_M. In turn, V_M is also equal to the area A₅ + A₆. Equating the two areas, we obtain A₇ = A₄ + A₅, which, as anticipated, divides the shaded regions in Fig. 3 B into equal parts. The fact that V_M 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 Δq_T = Q_max/N. Assuming that the value of Δq_T 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 V_M (which is independent of N) can be used to monitor the work function W_C[q] = −Δq_TV_M.

Despite the historical precedent for using the adjective “median” to describe V_M—arguments for which include not only the “equal areas” property but also the fact that when V = V_M, 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 V_M as a function of the normalized electrical capacitance f_q = Δq_T⁻¹d〈q〉/dV (Fig. 3 C):

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

Of the methods used to acquire f_q, 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 f_q 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 f_q in one sweep, is later demonstrated here through simulation. We note that f_q 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 V_M 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:

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

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 dV_M/dμ in BK channels, and knowing Δq_T, one can obtain a lower bound (= ^VΔ〈n〉_max) for the total number of calcium-binding sites Δn_T as an alternative to measuring the Hill coefficient.

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

The LPFs Z_V(−) and Z_V(+) in Eq. 32 correspond to sub-schemes of Scheme 2 derived from extreme potentials V₍₋₎ and V₍₊₎. Z_V(−) describes J and L particles as fixed in their resting positions, whereas K particles are allowed to bind Ca²⁺, and Z_V(+) is analogous except with J and L fixed in their activated positions. The LPFs contain CPFs Z_CR and Z_OA, as follows: Z_V(−) = Z_CR⁴ and Z_V(+) = L_V(+)(J_V(+)DZ_OA)⁴. Recalling that η_L = −kTlnL and η_J = −kTlnJ, we insert the above expressions for Z_V(−) and Z_V(+) into Eq. 32, yielding:

where we have eliminated the diverging second term in Eq. 32 by recognizing that, because Δq_T = Δq_L + 4Δq_J, the electrical components of η_L,V(+) and η_J,V(+) exactly cancel Δq_TV₍₊₎, leaving only the charge-neutral potentials η_L,V=0 = ΔG_L − Δn_Lμ and η_J,V=0 = ΔG_J − Δn_Jμ. In unusual cases where obligatory calcium binding is associated with the activation of L and/or J particles, one would observe sloping asymptotes (m = Δn_L + 4Δn_J) when plotting W_C[q] versus μ, but for standard models of BK, the global linkage plot has horizontal asymptotes.

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

The right side of Eq. 34 follows immediately from Eq. 33 by substituting the expression (1 + CEK)/(1 + K) for Z_OA/Z_CR, and recognizing that this is equal to the expectation value 〈CE〉_K. 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ΔW_E can be obtained by subtracting Eq. 34 from Eq. 23a.

It is reasonable to enquire, given the outcome of Eq. 34, whether a limited version of global analysis can be used to measure ΔW_E independently of ΔW_C. 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:

where V_M(C) and V_M(O) are mean activation potentials derived from the limiting curves 〈q_C〉 = kT(∂lnZ_L(−)/∂V)_μ and 〈q_O〉 = kT(∂lnZ_L(+)/∂V)_μ. Using Eq. 12, these can be expanded to:

In BK channels, it is possible to resolve 〈q_C〉 and 〈q_O〉 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 〈q_C〉 and 〈q_O〉 from measurements of G-V and Q-V plots, for example, through the relations 〈q_C〉 = 〈q〉 + kTdln(G_max − G)/dV and 〈q_O〉 = 〈q〉 + kTdlnG/dV. The expression kTdlnG/dV is known as the mean activation charge displacement 〈q_a〉, 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 〈q_a〉 (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 e_o), compatible with that of Δq_L. In contrast, the limiting slope observed in Shaker is substantially larger (12–13 e_o; see Islas and Sigworth, 1999)—equal in value to Δq_T 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 〈q_C〉 and 〈q_O〉. Starting from the identity W_H[g] = kTln(Z_L(+)/Z_L(−)), we derive dW_H[g]/dV = 〈q_O〉 − 〈q_C〉 (see also Conti, 1986), and recognizing that 〈q〉 = (1 − P_O)〈q_C〉 + P_O〈q_O〉, we obtain 〈q_C〉 = 〈q〉 − (G/G_max)dW_H[g]/dV and 〈q_O〉 = 〈q〉 + (1 − G/G_max)dW_H[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 ΔW_E for the simple reason that the J–K interaction does not directly involve the pore.

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 〈q_a〉 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 (m_V(−) = Δq_T) and positive (m_V(+) = Δq_L) 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):

Specifically, the voltage sensor and pore are both voltage sensitive (Δq_L, Δq_J ≠ 0), whereas the calcium-binding domain demonstrates no voltage sensitivity (Δq_K = 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.”

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.

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

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 Mg²⁺ (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 L→J_[G]→K_[F] is given by:

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:

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. 37a–d, 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.

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 (K₁ and K₂), for a total of 17 interacting gating particles. The particles L, J₂, K₁, and K₂ form an allosteric loop involving the core interactions C, D, and E, and an internal interaction F coupling K₁ and K₂. The J₁ particle does not participate in the loop but influences its neighboring J₂ 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→J₂→(J₁, K₂)→K₁ is given by:

from which we derive, as usual, the first-order equation Z = Z_C⁴ + LZ_O⁴. Expressions for the CPFs Z_C and Z_O are as follows:

The second-order CPFs Z_CR, Z_CA, Z_OR, and Z_OA are functions of K₁ and K₂, given by Z_CR = f₂(K₁, K₂) = (1 + K₁) + K₂(1 + K₁F), Z_CA = f₂(K₁E, K₂), Z_OR = f₂(K₁, K₂C), and Z_OA = f₂(K₁E, K₂C). The mean single-channel gating charge displacement 〈q〉 ranges from 0 to Δq_T = Δq_L + 4(Δq_J1 + Δq_J2) and is calculated from 〈q〉 = 〈l〉Δq_L + 〈j₁〉Δq_J1 + 〈j₂〉Δq_J2, where the voltage-dependent particle activation curves are derived as follows:

The average single-channel conductance is 〈g〉 = 〈l〉g_L.

Numerically computed equilibrium curves derived from Scheme 4 are shown in Figs. 6 (W_H[g] “Hill” analysis) and 7 (W_C[q] “capacitance” analysis) for a range of [Ca²⁺]_i values. Because Scheme 4 is a special case of the generalized Scheme 3, Eqs. 37a–d 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. 8–10, 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 ΔW_C + ΔW_E, 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 Q_max was obtained by integrating gating currents, slightly underestimated the true value of Δq_T = 2.62 e_o by −2.4%, a difference that is close to being statistically significant (P = 0.055).

The major challenge in generating W_H[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, V_eq, and g_leak. From the definition W_H[g] = kTlnG − kTln(G_max − G), we recognize that the kTlnG defines the V₍₋₎ asymptote of W_H[g] and is logarithmically sensitive to the floor, whereas kTln(G_max − G) 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 = Ng_LP_O and variance Var_G = Ng_LP_O(1 − P_O) through the relation S/N = G/Var_G^1/2, which works out to (NG/(G_max − G))^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 P_O < 1/N. Therefore, with a cumulative N of 10⁵ (= 10² traces × 10³ channels), noise effects became evident for P_O ≤ 10⁻⁵, as seen in Fig. 9. It is remarkable that useful data could be derived for P_O smaller than the predicted limit, down to P_O = 10⁻⁸. The corresponding S/N ratio for the ceiling variable G_max − G is (N(G_max − G)/G)^1/2. Thus, for S/N ≈ 1, and cumulative n = 10⁵, noise is increased when 1 − P_O = 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 P_O, it was least affected by floor and ceiling effects and was used to fit the slope Δq_L. The average fitted value of Δq_L differed only slightly (+1.6%) from the given value of 0.3 e_o. To determine interaction energies ΔW_C and ΔW_D, 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, V_eq, and g_leak) so that they visually paralleled the (V₍₋₎, μ₍₊₎) asymptote. The value of ΔW_C (averaged from measurements taken at V₍₋₎ and V₍₊₎) was 1.3% less than the given value of −53 meV. Measuring ΔW_D 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” Ca²⁺ plot was excluded from the analysis.

In W_C[q] experiments (Fig. 10), the average value of the outcome parameter ΔW_C + ΔW_E 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 W_C[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 [Ca²⁺]_i. It was quantified by simulating the model in the absence of Nyquist noise, and the ^VΔV_M 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 Δq_T measurements, accounted for most of the error in ΔW_C + ΔW_E.

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 (J₁, J₂, K₁, and K₂) in Scheme 4 are shown in Fig. 11. The work functions for the Ca²⁺-sensitive K particles are plotted on μ and log[Ca²⁺]_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 J₁ providing the best resolution because of its small intrinsic charge and relatively large interaction energy with J₂. From an experimental standpoint, the Hill plot for J₂ 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 K₁ and K₂ (ΔW_F > 0).

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 B₁, B₂, B₃, and B₄, 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 ΔG_P = ΔG_L/4 was additionally offset by the amount −ΔW_B = −(ΔW_B1 + ΔW_B2 + ΔW_B3 + ΔW_B4)/4 to maintain the original free energy change ΔG_L 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 −ΔW_B and −2ΔW_B; Fig. 12 B) and brief occupancy times. These attributes of the interacting tetrameric pore are found in the Scheme 5 partition function:

where Z_C and Z_O remain unchanged from Scheme 4 (Eq. 38). It is readily apparent that for large values of B = exp(−ΔW_B/kT), the partition function reverts back to Scheme 4 with its binary pore equilibrium constant L = P⁴. 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 ΔW_B = −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: B_p = B^1/4, 〈g〉 = 〈p〉g_L/4, 〈p〉 = ∂lnZ/∂lnP, and η_P = η_L/4. The results of global capacitance (W_C[q]) analysis were practically unchanged from Scheme 4 because interactions between voltage- and Ca²⁺-sensitive particles are unaffected by B; therefore, W_C[q] plots are not shown. However, dramatic decreases in the limiting slopes of W_H[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 W_H[g] were equal to:

The rise ^VΔ(W_H[g] − η_P)_μ(±) in asymptotes taken from the difference of Eqs. 41a and 41b equals −(ΔW_D + 2ΔW_B). The value of −2ΔW_B 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 ΔW_B 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 B_p and subconductance levels g_p = x_pg_L. The loss of uniform incremental values in these parameters influences the measurement of self-interaction energy, as demonstrated by the new linkage relation:

The second term on the right of Eq. 42 is an artifact of the nonuniform spacing of subconductance levels and vanishes if x₁ = 0.25 and x₃ = 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(B³/B₁B₂), the difference of which is kTln(B₃B₄) = −(ΔW_B3 + ΔW_B4). These expressions stem from the constraints imposed on self-interaction energies to achieve a “balanced” equilibrium between closed and open states.

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 Ca²⁺-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 e_o. Another mechanism for dual sensitivity is the obligatory binding of Ca²⁺ during pore opening. As unlikely as this seems, Piskorowski and Aldrich (2002) have reported Ca²⁺ 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 Ca²⁺ upon opening was performed by specifying the particle function μ_L = ΔG_L + −Δq_LV − Δn_Lμ, where Δn_L was chosen to be 0.4.³ This led to sloping asymptotes in the W_C[q] versus μ plot (Fig. 13 A). The energy Hill plot W_H[g] demonstrated asymptotic slopes for both V and μ axes (Fig. 13 B).

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, Mg²⁺ 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 Mg²⁺ displacement could affect the interaction energy ΔW_E. 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 ΔW_G = ΔG_G − Δq_GV. A voltage-sensitive G = exp(−ΔG_G/kT) would contribute an amount (∂lnZ/∂lnG)Δq_G of gating charge to the Q-V curve. The effects on the Q-V and Hill plots from a small charge (Δq_G = 0.2 e_o) are shown in Fig. 14 (A and B). Because interactions between voltage- and Ca²⁺-sensitive components of the channel do not depend on internal energies within the voltage sensor, global analysis (Fig. 14 A) was unaffected except that Δq_T grew larger by the amount 4Δq_G. On the other hand, the V asymptotes in the local Hill plot for J₂ diverged because of the addition of Δq_G to the slope of the upper asymptote (Fig. 14 B). The interaction energy −(ΔW_G + ΔW_D) 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 Δq_E = 0.2 e_o was assigned to the allosteric factor E, which, unlike G, couples particles of different modalities (J₂ and K₁). Here, a more dramatic change to the Q-V curve occurred (Fig. 13 C). At very low [Ca²⁺]_i, where E ≈ 0, the channel behaved as in the original Scheme 4 with a well-defined V_Mμ(−), despite the theoretical value of V_M being shifted immeasurably far to the right on the V axis. With increasing [Ca²⁺]_i, The E charge component (∂lnZ/∂lnE)Δq_E entered the Q-V plot from the right side, with V_M eventually reaching a well-defined limiting value specified by V_Mμ(+). The global energy of interaction −4(ΔW_C + ΔW_E) can be deduced from the area A₂ − A₁ shown in Fig. 14 C. From an experimental standpoint, obtaining the correct result would require that one recognize that Δq_T, the total charge movement per channel, has become a function of [Ca²⁺]_i as a result of linkage between J and K. Complex behavior was also demonstrated in the W_H[j2] plot (Fig. 14 D), which contained numerous inflection points, generating “false” and ultimately diverging asymptotes across a wide range in voltage.

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 W_C[q], a global parameter derived from the Q-V plot; and the Hill energy W_H, 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 W_H[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.

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.

The major source of error in W_C[q] measured in voltage-ramp simulations was systematic, caused by a Ca²⁺-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 V_M 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 >10⁶ 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 V_M “creep”) and an environmental dependence in the value of Q_max, 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 (= mC_patch) 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⁻⁴ e_o/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 ΔC_J = ∂Δq_J/∂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 ∂Δq_J (assuming ΔC_J is voltage independent) would obtain a second-order voltage term: η_J = ΔG_J − Δq_oJV − ΔC_J(V − V_o)V. This could theoretically have an effect on baselines and limits of integration when analyzing ramp-generated C_q to measure V_M. 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 W_C[q] requires knowing the value of total gating charge, the so called “charge-per-channel” Δq_T = Q_max/N. Although a model-dependent estimate of Δq_T = 2.62 e_o has been determined for BK (Horrigan and Aldrich, 2002), an independent determination of Δq_T 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).

The critical step in W_H[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 P_O < 1/N for the floor asymptote, and 1 − P_O = 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 P_O (“ceiling”) asymptote of W_H[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⁻⁷ (see Fig. 9), requiring at the very least a 24-bit A/D converter (10⁷ ≈ 2²⁴). An elegant solution to this problem would be to construct a “Hill” circuit, which generates analogue signals for G and G_max − G by adjusting the driving force V − V_eq and baseline I_base of the ionic current, and then applies a logarithmic converter before digitizing. Using a repeating ramp protocol, manual or semi-automated adjustment of inputs G_max, V_eq, and I_base 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 P_O and 1 − P_O. 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).

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 Ca²⁺ dependence of P_O at very negative potentials, the authors computed a variable they called R_O that is equivalent to Z_OR/Z_CR, and by plotting logR_O versus [Ca²⁺] obtained log〈C〉. They also obtained the intrinsic pore charge Δq_L by fitting the negative asymptote of log(P_O). They recognized the importance of using log(P_O), 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 P_O (using an all-points histogram of single-channel currents to measure NP_O and dividing by N ≈ G_max/g_L) 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 W_H[g] versus V. Instead, an alternative method for estimating D using the product of Δq_J and the midpoint voltages of 〈q_C〉 and 〈q_O〉 was proposed specifically for the HA model (Scheme 2). The proof for the more general Scheme 3 goes as follows: recognizing that 〈q_O〉 − 〈q_C〉 = dW_H[L]/dV, we integrate both sides with respect to V to obtain ∫〈q_O〉dV − ∫〈q_C〉dV = ^VΔW_H[g]. Separating 〈q_O〉 into its baseline (Δq_L) and variable (〈q_O〉 − Δq_L) components, we subtract ^VΔμ_L = Δq_L(V₍₊₎ − V₍₋₎) from both sides of the equation and integrate the remaining expression on the left side by parts, obtaining −Δq_J(V_M(O) − V_M(C)) = ^VΔ(W_H[g] + η_L). V_M(C) and V_M(O) are the first moments of normalized capacities derived from 〈q_C〉 and 〈q_O〉 − Δq_L, respectively. Recalling that ^VΔ(W_H[g] + η_L) = −4kTln(D), and defining ^LΔV_M = V_M(O) − V_M(C), we obtain the desired result: Δq_J ^LΔV_M = 4ΔW_D (see Fig. 15). Assuming one can measure both 〈q_C〉 and 〈q_O〉—for example, by integrating the fast component of gating charge after prepulsing to either V₍₋₎ or V₍₊₎—then ΔW_D could be calculated without resorting to logarithmic (Hill) analysis of the G-V curve. The value of Δq_J, which is necessary to complete the calculation, can theoretically be obtained by subtracting Δq_L from Δq_T, both of which can be independently measured. The authors acquired 〈q_O〉 indirectly through the relation 〈q_O〉 = 〈q_a〉 + 〈q〉, where 〈q_a〉 = 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(P_O)-V curves.

The allosteric factor E was estimated in similar fashion by measuring shifts in 〈q_C〉 and 〈q_O〉 in response to saturating changes in [Ca²⁺]_i. Because the influence of pore activation has been removed from both 〈q_C〉 and 〈q_O〉, 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 ΔW_G 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 P_O was very low (∼10⁻⁷), 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 NP_O histogram to generate positive voltage (“ceiling”) asymptotes in the conductance Hill plot.

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/(M_max − M)), 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 + Kⁿ, 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 = Kⁿ and whose activation curve is therefore Kⁿ/(1 + Kⁿ), 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 n_H 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 Δq_B = n_HΔ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 V_M and the expression W_C[q] = −Δq_TV_M, 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 V_1/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, V_1/2 is often used as a stand-in for the more correct V_M. In a somewhat confusing twist of nomenclature, the V_1/2 (the voltage for which the normalized Q-V is 0.5) is equal to the median value of the capacitance distribution function f_q (see accompanying article by Chowdhury and Chanda in this issue), whereas V_M, 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, V_1/2 has no special thermodynamic significance (in the sense of describing an energy or displacement), but it may coincide with V_M for symmetrical Q-V plots.

Similarly, the product Δq_BV_1/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〉 = Δq_T − 〈q_a〉, where 〈q_a〉 = kTdlnG/dV is the mean activation charge displacement. Therefore, gating current measurements are not required to measure W_C[q], because even Δq_T can be obtained from the limiting value of activation charge displacement: Δq_T = 〈q_a〉_V(−). However, the logarithmic dependence of 〈q_a〉 on G implies that the extreme edges or “tails” of the G-V are important in evaluating W_C[q], which for a late-activating Kv channel is related to G through −kT∫V(d²lnG/dV²)dV. This contrasts with the traditional Boltzmann method of analysis, where Δq_BV_1/2 obtained by Boltzmann fitting is numerically equal to −kTln(B/(1 − B))_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 Δq_BΔV_1/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 W_C[g] = Δq_L∫VdP_O = Δq_LV_M[g], where V_M[g] = ∫Vf_gdV and f_g = dP_O/dV. The variable V_M[g] could be referred to as the “mean” activation voltage of conduction, and f_g 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 Δq_L^μΔV_M[g], as a kind of “global analysis” of the local K–L interaction, is equal to −4ΔW_C. This attempt at deriving the strength of K–L linkage (or any perturbation affecting pore opening, say from a mutation) by measuring shifts in V_M[g] is flawed, because, unlike Q = NkT(∂lnZ/∂lnV), the conductance G = (∂lnZ/∂lnL)G_max 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 W_C[g] by parts to obtain Δq_L∫P_OdV − Δq_LV₍₊₎, and then using the chain rule on the first term, as in:

are problematic, particularly as the second term diverges for Δn_L = 0. Recognizing that P_O = −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 −ΔW_E (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 Δq_BV_1/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.

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 = ΔH_K − ΔS_KT, 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 T_o = ΔH_K/ΔS_K. The steepness of the curve is determined by ΔS_K. Rearranging K to read

we can estimate the proportional increase in the density of states, R_Ω = exp(ΔS_K/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 ∼10⁵⁵, which is equivalent to an entropy change of ΔS_K = 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 ΔH_K = ΔS_KT_o ≈ −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 ΔW_C. The enthalpic component of ΔW_C can be measured from the rise of the sigmoidal component of the conductance Hill energy W_H[g] plotted versus T. If there is also an entropic component ΔS_C, it would show as a difference in the slopes of the positive and negative asymptotes equal to ΔS_C.

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 ΔC_p experienced by a protein upon denaturing is typically about five times the value of ΔS_Ko evaluated at T_o (Privalov and Khechinashvili, 1974), which would generate a measurable nonlinearity in the particle potential η_K of a T-sensor. The value of ΔC_p appears to be fairly constant (Baldwin, 1986), leading to η_K = (ΔC_p − ΔS_Ko)(T − T_o) − ΔC_pTln(T/T_o), which is derived by integrating the thermodynamic identities C_p = (∂H/∂T)_P = T(∂S/∂T)_P around T_o = ΔH_Ko/ΔS_Ko, 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 W_H[g] related to the possibility that η_K crosses zero twice, theoretically generating two melting temperatures separated approximately by 2ΔT = 2ΔH_Ko/ΔC_p, 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 T_o, 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.

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 >10³ 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 〈q_a〉 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)⁴ + L(1 + KC)⁴, and for Scheme 2, it is Z = [(1 + K) + J(1 + KE)]⁴ + L[(1 + KC) + JD(1 + KCE)]⁴. 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 Ca²⁺ sensors K₁, K₂, K₃… can be represented by K_[F], where F = {F₁, F₂, F₃…} 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])₍₊₎ = (K₁K₂K₃…F₁F₂F₃…). 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.