A thermodynamic approach to studying allosterically regulated ion channels such as the large-conductance voltage- and Ca^{2+}-dependent (BK) channel is presented, drawing from principles originally introduced to describe linkage phenomena in hemoglobin. In this paper, linkage between a principal channel component and secondary elements is derived from a four-state thermodynamic cycle. One set of parallel legs in the cycle describes the “work function,” or the free energy required to activate the principal component. The second are “lever operations” activating linked elements. The experimental embodiment of this linkage cycle is a plot of work function versus secondary force, whose asymptotes are a function of the parameters (displacements and interaction energies) of an allosteric network. Two essential work functions play a role in evaluating data from voltage-clamp experiments. The first is the conductance Hill energy *W _{H}*

_{[g]}, which is a “local” work function for pore activation, and is defined as

*kT*times the Hill transform of the conductance (G-V) curve. The second is the electrical capacitance energy

*W*

_{C}_{[q]}, representing “global” gating charge displacement, and is equal to the product of total gating charge per channel times the first moment (

*V*) of normalized capacitance (slope of Q-V curve). Plots of

_{M}*W*

_{H}_{[g]}and

*W*

_{C}_{[q]}versus voltage and Ca

^{2+}potential can be used to measure thermodynamic parameters in a model-independent fashion of the core gating constituents (pore, voltage-sensor, and Ca

^{2+}-binding domain) of BK channel. The method is easily generalized for use in studying other allosterically regulated ion channels. The feasibility of performing linkage analysis from patch-clamp data were explored by simulating gating and ionic currents of a 17-particle model BK channel in response to a slow voltage ramp, which yielded interaction energies deviating from their given values in the range of 1.3 to 7.2%.

## INTRODUCTION

Ion channels are pore-forming proteins that regulate (gate) the flux of selected ions (K^{+}, Na^{+}, Ca^{2+}, 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^{2+}-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^{2+} 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^{2+}-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^{2+}]_{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^{2+}-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*. 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”

_{YX}*η*, 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*and not the total energy of activating both components (=

_{AB}*W*+

_{A}*W*+

_{B}*W*). This is helpful because, in many instances, it is not possible to measure

_{AB}*W*, 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

_{B}*F*such as voltage or calcium chemical potential as secondary agents, as these can be smoothly and quantitatively controlled in the laboratory.

_{B}Using a combination of lever pulls, the energy of interaction *W _{AB}* between components A and B can be obtained by twice measuring

*W*(i.e., pulling the A lever), once with component B in the resting state (

_{A}*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*, the logic being that if A and B are allosterically linked, activating B will assist (or retard) activation A by the amount

_{AB}*W*. This form of linkage analysis can be concisely summarized by the expression

_{AB}*W*=

_{AB}*Δ*

^{B}*W*=

_{A}*Δ*

^{B}*ΔΦ, where the “lever operator” Δ() is defined (using component A as an a example) by*

^{A}*Δ() = ()*

^{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*, whereas the other leg is defined by the lever operation

_{A}*Δ (see Fig. 2 B). The reciprocity principle demands that the alternative cycle, which uses*

^{B}*W*and

_{B}*Δ as its legs, should yield the same interaction energy; in other words,*

^{A}*W*=

_{AB}*W*. 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

_{BA}*W*(in which the activation state of a third component C influences

_{ABC}*W*), are not covered here but represent a straightforward extension of the theory.

_{AB}The allosteric relationship between components A and B is established by plotting *W _{A}* against

*F*, referred to here as a linkage plot. The typical appearance of a linkage plot is shown in Fig. 1. The work function

_{B}*W*rises steeply in the region of B activation, with a displacement equal to

_{A}*W*. If the principal component A is itself not sensitive to

_{AB}*F*, the plot plateaus at the two extremes. Otherwise, the extreme behavior is characterized by sloping asymptotic lines whose derivative

_{B}*m*is the generalized displacement of A linked to

*F*(for example, if

_{B}*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*is derived from the G-V curve through

_{L}*W*= −

_{L}*kT*ln[

*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*will be referred to as the “conductance Hill energy”

_{L}*W*

_{H}_{[g]}. Defining

*W*

_{H}_{[g]}as the negative value of −

*W*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,

_{L}*W*and

_{J}*W*, which are the work functions for voltage and Ca

_{K}^{2+}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*is derived by integrating

_{q}*Vdq*over the saturating range of

*q*given by [0, Δ

*q*]. Changing the integrand from

_{T}*q*to

*V*yields

*W*= Δ

_{q}*q*, where

_{T}V_{M}*V*= ∫

_{M}*Vf*is the mean value of

_{q}dV*V*with respect to the capacitance distribution

*f*= (

_{q}*dq*/

*dV*)/Δ

*q*. Because the work function for

_{T}*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*, 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,

_{q}*W*= Δ

_{n}*n*describes the reversible energy of ligand binding, where

_{T}μ_{M}*μ*is the mean value of the chemical potential

_{M}*μ*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

^{2+}-sensing elements in the BK channel is equal to the expression

*Δ*

^{μ}*W*= Δ

_{q}*q*Δ

_{T }^{μ}*V*, whose evaluation in practice involves measuring the shift in

_{M}*V*under limiting conditions of “zero” and “saturating” calcium concentration, denoted by

_{M}*μ*

_{(−)}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*can be obtained. Wyman’s reciprocity principle assures us that

_{M}*Δ*

^{V}*W*=

_{n}*Δ*

^{μ}*W*, or more succinctly, that

_{q}*W*=

_{qn}*W*. In other words, measuring the mean ligand activity

_{nq}*μ*as a function of

_{M}*V*could in principle be used to confirm the value of

*W*obtained from the

_{qn}*μ*dependence of

*V*, provided one is able to measure Ca

_{M}^{2+}-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^{2+}-binding site located a small electrical distance *d* inside the membrane potential. If the resulting charge movement 2*d* 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

^{2+}bound per site. However, in most models of BK, there is the assumption of strict modality separation between voltage and Ca

^{2+}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.

## MATERIALS AND METHODS

### Data simulations

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

*α*. During constant-voltage segments of the simulation, or if the particle charge Δ

_{i}*q*= 0,

_{i}*d*was obtained from the formula −ln

_{i}*r*/

_{n}*α*, where

_{i}*r*is a uniform random number drawn between 0 and 1 obtained with the long-period (>2 × 10

_{n}^{18}) random number generator of L’Ecuyer with Bays-Durham shuffle (Press et al., 1992). During the ramp phase, the

*d*of charged particles was obtained from the formula (2

_{i}*kT*/

*m*Δ

*q*)ln[1−(

_{i}*m*Δ

*q*ln

_{i}*r*)/(2

_{n}*kTα*

_{oi})] (Sigg and Bezanilla, 1997), which reflects the nonstationary value of the forward and backward particle rate constants

*α*and

_{i}*β*, which are of the form

_{i}*α*

_{oi}exp(Δ

*q*/2

_{i}V*kT*) and

*β*

_{oi}exp(−Δ

*q*/2

_{i}V*kT*), respectively. The method for determining

*α*and

_{i}*β*are given in the next section. The particle with the shortest

_{i}*d*was allowed to transition, after which the simulation time was advanced by

_{i}*d*. If between

_{i}*t*and

*t*+

*d*there was a change in pulse protocol occurring at

_{i}*t*

_{1}, the simulation time reverted to

*t*

_{1}, 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*) using a Gaussian filter and discretized (sampling frequency

_{c}*fs*= 5

*f*) in such a way that the discrete integral of the response equaled the gating charge displacement Δ

_{c}*q*(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/

_{i}*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*and driving force (

_{L}*V*−

*V*). The value of the equilibrium potential

_{eq}*V*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 (

_{eq}*I*) was equal to the Nyquist value 4

_{noise}*kTBg*, where

*B*is the low-pass bandwidth of a Gaussian filter (

*B*≈ 1.064

*f*; see Crouzy and Sigworth, 1993), and the conductance

_{c}*g*was made up of contributions from a constant “leak” conductance

*g*, assumed to be 50.0 ± 2.5 pS (≈ 20 Gohm

_{leak}^{−1}), and a fluctuating term

*N*, where

_{O}g_{L}*N*was the number of open channels at a given time. Each simulated current trace (gating and ionic) contained the following terms:

_{O}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*, and

_{leak}*C*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

_{patch}*V*= −300 mV, ending with a tail (OFF) pulse to 0 mV.

_{r}### Kinetics of the 17-particle model

Although the emphasis of this paper is on equilibrium methods, kinetic parameters for the 17-particle model were required to perform ramp simulations. The chosen scheme was consistent with known macro-kinetics of the channel while preserving equilibrium properties. Time constants of decay *τ* for both ionic and gating currents in BK were obtained from Horrigan and Aldrich (2002). In their study, maximum value *τ _{max}* of the voltage-sensor decay time, estimated from the

*τ*versus

*V*plot of the rapid component of gating current, was ∼5 ms, whereas ionic current kinetics, reflecting pore activation, demonstrated a maximum

*τ*in the decay time of ∼0.05 ms. Both voltage-sensor particles in Scheme 3 were assigned the same decay time. Ca

_{max}^{2+}binding was assumed to have the same

*τ*as voltage-sensor activation. The condition of maximum

_{max}*τ*= (

*α*+

*β*)

^{−1}was obtained when forward (

*α*) and backward (

*β*) particle activation rates were both equal to (2

*τ*)

_{max}^{−1}. To satisfy equilibrium and kinetic requirements, particle rate constant were computed according to

*α*= (2

_{i}*τ*)

_{max}^{−1}(

*Z*

_{i}_{(+)}/

*Z*

_{i}_{(−)})

^{1/2}and

*β*= (2

_{i}*τ*)

_{max}^{−1}(

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

### Data analysis

Simulated experiments consisted of 100 traces each of ionic and gating currents from *n* = 1,000 ± 50 channels, taken at [Ca^{2+}]_{i} = 10^{−3} µM and 10^{3} µM. Estimates for *g _{leak}* and

*C*were obtained by analyzing short segments (2–40 ms) of current before (

_{patch}*I*

_{(−)}) and after (

*I*

_{(+)}) the beginning of the test pulse (

*t*

_{1}). 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*and

_{leak}*C*were simultaneously obtained by fitting

_{patch}*I*to the following function:

*I*(

*t*<

*t*

_{1}) =

*I*

_{(−)};

*I*(

*t*≥

*t*

_{1}) =

*I*

_{(−)}+

*mC*+

_{patch}*g*(

_{leak}*V*(

*t*) −

*V*). Fits were obtained using the solver application in Excel (Microsoft). The values of

_{r}*g*and

_{leak}*C*, in addition to

_{patch}*V*(in ionic current experiments), were further refined through manual adjustment and fitting of asymptotes to achieve satisfactory shapes for both the Q-V curve = ∫(

_{eq}*Ni*/

_{g}*m*)

*dV*and the conductance Hill plot

*W*

_{H}_{[g]}=

*kT*ln(

*G*/(

*G*−

_{max}*G*), where

*G*=

*I*/(

*V*−

*V*). To reduce the deleterious effects of Nyquist current noise on Q-V and

_{eq}*W*

_{H}_{[g]}plots, post-filtering of currents was performed after initial correction for

*g*and

_{leak}*C*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

_{patch}*V*was estimated from the following integral (see Eq. 29):

_{M}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

^{2+}]

_{i}= 10

^{−3}and 10

^{3}µM. The pore gating charge Δ

*q*was obtained from the slope of the lower (negative)

_{L}*V*asymptote in the high Ca

^{2+}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*

^{2}/

*N*(Sigworth, 1980) for both test and tail currents, with

*N*averaged over the two estimates. The total gating charge Δ

*q*was obtained from the integrated positive gating current Δ

_{T}*Q*= ∫

*Ni*divided by the measured

_{g}dV*n*value.

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

### Online supplemental material

The partition function and linkage relationships for Scheme 3 (see Eqs. 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 at http://www.jgp.org/cgi/content/full/jgp.201210859/DC1.

## RESULTS

### Thermodynamics

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

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

*μ*are the chemical potentials of buffered and fixed-numbered species present in numbers

_{f}*N*and

_{b}*N*. 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.

_{f}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*− Σ

*(notation adapted from Alberty, 2002). Integrating Eq. 1 at constant values of the intensive variables, we obtain G[*

_{b}μ_{b}N_{b}*V*,{

*μ*}] = Σ

_{b}*, 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[*

_{f}μ_{f}N_{f}*V*,{

*μ*}] =

_{b}*Nμ*. Following Wyman (1975), we equate the system energy Φ with

_{BK}*μ*. Writing the Gibbs–Duhem equation for Eq. 1 as −

_{BK}*SdT*+

*VdP*−

*QdV*− Σ

*−*

_{b}N_{b}dμ_{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^{2+} 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^{2+} chemical potential *μ* = *μ*^{ο} + *kT*ln*a _{Ca}* is a function of internal calcium activity

*a*(in practice, assumed to be [Ca

_{Ca}^{2+}]

_{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).^{1} 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**,

**, and**

*k**l*are the activation numbers for protomers

**J**,

**K**, and L, ranging from 0 to a maximum of 4 (

*n*,

_{J}*n*) or 1 (

_{K}*n*). The variables pertaining to the voltage sensor (=

_{L}**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}*η*

_{J}_{1}+ ∂Φ/∂

*η*

_{J}_{2}+…), and a product with composite allosteric factors (for example

**=**

*J**J*

_{1}

*J*

_{2}…).

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”

*η*as the change in Φ

_{X}*required to activate the X protomer (*

_{X}*η*≡

_{X}*ΔΦ*

^{X}*), and evaluating*

_{X}*η*at constant

_{X}*T*,

*P*,

*V*(voltage), and

*μ*, we obtain:

where Δ*q _{X}* and Δ

*n*are charge and Ca

_{X}^{2+}displacements of activation, and Δ

*G*= Δ(

_{X}*U/N)*− Δ

_{X}*s*+ Δ

_{X}T*v*is the Gibbs energy of particle activation (not to be confused with the macroscopic conductance

_{X}P*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*, because the pore is assumed to be purely voltage dependent.

_{L}VAn “accessible” allosteric system is defined as one for which the weights Θ_{j}_{,k,l} in Eq. 5 are nonzero for any combination of ** j**,

**, and**

*k**l*. The partition function

*Z*can be converted into polynomial form:

where the particle equilibrium constants ** J**,

**, and**

*K**L*are related to their respective particle potentials through “local bridge” equations (compare with Eq. 4)

*η*= −

_{J}*kT*ln

**,**

*J**η*= −

_{K}*kT*ln

**, and**

*K**η*= −

_{L}*kT*ln

*L*. 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**〉 = ∂Φ/∂

*η*= ∂ln

_{J}*Z*/∂ln

**is the mean number of activated**

*J***J**particles, and the other activation numbers are similarly given by 〈

**〉 = ∂Φ/∂**

*k**η*= ∂ln

_{K}*Z*/∂ln

**and 〈**

*K**l*〉 = ∂Φ/∂

*η*= ∂ln

_{L}*Z*/∂ln

*L*. 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**〉, 〈

**〉, and 〈**

*k**l*〉 are easily derived from the polynomial form of the partition function, Eq. 7.

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

The words “protomer” and “particle,” referring to the basic building blocks of allosteric networks, have been used almost interchangeably so far in this paper. However, it is useful to propose a distinction between the two terms to clarify the role of complex domains with multiple activation states. The word “particle” will hereafter refer to a single activated protomer state. In the case of a two-state protomer, such as the pore with partition function *Z* = 1 + *L*, the terms “protomer” and “particle” are synonymous. However, a composite protomeric unit possesses more than one activated state. This may be because it undergoes successive transitions during activation, as is thought to be true for the S4 segment in Kv channels (Schoppa and Sigworth, 1998; Gamal El-Din et al., 2010; Henrion, et al., 2012), in which case the partition function is written *Z* = 1 + *J*_{1} + *J*_{2} +…; or it is composed of multiple interacting two-state particles, in which case the partition function is (for two particles) *Z* = (1 + *J*_{1}) + *J*_{2}(1 + *J*_{1}*G*), with *G* as a coupling factor. In both cases, the particles are connected to the equilibrium constants *J*_{1}, *J*_{2}, 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*_{2}, 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*_{1}〉 = ∂ln*Z*/∂ln*J*_{1}, 〈*j*_{2}〉 = ∂ln*Z*/∂ln*J*_{2}, 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*^{4} from which the classical Hill equation (Hill, 1910) can be derived as *H* = (1/4)∂ln*Z*/∂ln*P* = *P*^{4}/(1 + *P*^{4}). The two-state partition function *Z* = 1 + *L* generates the Boltzmann equation, *B* = ∂ln*Z*/∂ln*L* = *L*/(1 + *L*), from which we conclude that in the concerted activating pore, the equilibrium constant is *L* = *P*^{4}.

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

### Contracted partition functions (CPFs)

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

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*and

_{C}*Z*, 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.

_{O}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^{2+}-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*= (1 +

_{O}*KC*). The ratio of the two CPFs yields the expectation value of C derived from the distribution of K particle states:

The Ca^{2+} dependence in Eq. 13 arises from the equilibrium constant *K* = exp(−*η _{K}*/

*kT*) = exp[(−Δ

*G*+ Δ

_{K}*n*)/

_{K}μ*kT*]. Recalling that [Ca

^{2+}]

_{i}= exp(

*μ*/

*kT*), this can be rearranged to read

*K*= ([Ca

^{2+}]

_{i}/

*K*

_{d})

^{ΔnK}, where

*K*

_{d}= exp(Δ

*G*/Δ

_{K}*n*) is the intrinsic disassociation constant of the K-binding site. If a single ion is bound during activation, then Δ

_{K}kT*n*= 1, and we obtain the usual expression

_{K}*K*= [Ca

^{2+}]

_{i}/

*K*

_{d}. Because

*K*grows monotonically with increasing

*μ*, it is correct to substitute 〈

*C*〉

*for 〈*

_{μ}*C*〉

*on the right side of Eq. 13. The value of 〈*

_{K}*C*〉

*changes smoothly from 1 to*

_{μ}*C*when varying

*μ*over the saturating range [

*μ*

_{(−)},

*μ*

_{(+)}]. Switching to energy units, we conclude that

*kT*ln〈

*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 conductance Hill energy W_{H[g]}

The Hill energy is defined as the work function of a single particle. It was stated in the Introduction that the pore-specific Hill energy *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}*kT*ln

*Z*

_{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*of channel opening as

_{O}*G*/

*G*=

_{max}*Z*

_{L}_{(+)}/

*Z*. Combining the above equations, and recalling that the conductance Hill energy was defined as

*W*

_{H}_{[g]}= −

*W*, we obtain, as advertised, the following Hill equation:

_{L}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]}+

*η*varies in sigmoid fashion from 0 to −4Δ

_{L}*W*in the range [

_{C}*μ*

_{(−)},

*μ*

_{(+)}]. 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^{2+} dependent, its particle potential *η _{L}* = Δ

*G*− Δ

_{L}*q*is independent of

_{L}V*μ*; i.e.,

*Δ*

^{μ}*η*= 0. Under these circumstances,

_{L}*η*can be left out of Eq. 16.

_{L}Plotting *W _{H}*

_{[g]}versus

*V*instead of

*μ*yields a straight line given by

*W*

_{H}_{[g]}= −

*η*. This is the Hill plot of the pore particle in the absence of interactions with other voltage sensors. Under conditions of “zero” Ca

_{L}^{2+}, symbolized by

*μ*

_{(−)}, the line has slope

*m*= Δ

*q*and

_{L}*V*-intercept

*V*

_{L,μ}_{(−)}= Δ

*G*/Δ

_{L}*q*, consistent with

_{L}*η*= −Δ

_{L}*q*(

_{L}*V*−

*V*

_{L,μ}_{(−)}). If we increase Ca

^{2+}to saturating levels, then

*m*is unchanged but the

*V*-intercept

*V*

_{L,μ}_{(+)}= (Δ

*G*+ 4Δ

_{L}*W*)/Δ

_{C}*q*is left-shifted in the setting of positive allosteric interaction (Δ

_{L}*W*< 0). This single-particle behavior is the source of sloping asymptotes in Hill plots.

_{C}### The ligand-binding Hill plot for Scheme 1

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

*Nb*)? 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*〉 = ∂ln

*Z*/∂ln

*K*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:

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*(however, notice the absence of the factor 4 compared with Eq. 16). Because the particle potential

_{C}*η*= Δ

_{K}*G*− Δ

_{K}*n*is a linear function of

_{K}μ*μ*, the plot of

*W*

_{H}_{[b]}versus

*μ*demonstrates sloping asymptotes (

*m*= Δ

*n*), as illustrated in Fig. 1 (substituting

_{K}*W*

_{H}_{[b]}for

*W*, and

_{A}*μ*for

*F*). Assuming Δ

_{B}*n*= 1, the expressions for the

_{K}*μ*intercepts of the two asymptotic lines are

*μ*

_{K,μ}_{(−)}= Δ

*G*and

_{K}*μ*

_{K,μ}_{(+)}= Δ

*G*+ Δ

_{K}*W*, from which the disassociation constant

_{C}*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*< 0).

_{C}### Interacting K particles in Scheme 1

Thus far, we have assumed that each K particle is unaware of the state of the other three K particles except as inferred from the activation state of the pore. In reality, there is extensive intersubunit contact between Ca^{2+}-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*

_{0}Z_{(0)},

*B*

_{0}B_{1}Z_{(1)},

*B*

_{0}B_{1}B_{2}Z_{(2)}, 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

_{0}*kT*ln

*B*and

_{1}*kT*ln

*B*to Eqs. 21a and 21b, respectively. Therefore, in evaluating

_{4}*Δ*

^{μ}*W*

_{H}_{[b]}for Scheme 1, we obtain—in addition to the K–L interaction energy −Δ

*W*and the divergent term −

_{C}*Δ*

^{μ}*η*= Δ

_{K}*n*(

_{K}*μ*

_{(+)}−

*μ*

_{(−)})—a third, self-interaction term

*kT*ln(

*B*/

_{4}*B*), representing the difference in the work needed to activate the first and last K particles.

_{1}### Summary of energy Hill plot analysis

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

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

### Scheme 2 (the HA model)

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

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*, derived from the parenthesized contents in the linkage diagram, are now functions of

_{O}*J*and

*K*(Table 4). A single polynomial function

*f*

_{1}describes both CPFs:

*Z*=

_{C}*f*

_{1}(

*J*,

*K*,

*E*) and

*Z*=

_{O}*f*

_{1}(

*JD*,

*KC*,

*E*). Depending on whether

*J*or

*K*is varied, the CPF ratio

*Z*/

_{O}*Z*can be interpreted as an expectation value for either

_{C}*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}*Δ↔

*Δ and*

^{V}*Δ↔*

^{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*and Δ

_{D}*W*are nonzero. Otherwise, J could be viewed as merely an accessory component of either the pore or the Ca

_{E}^{2+}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*(see Eq. 22a and Table 4). Therefore,

_{O}*E*must be determined from the Q-V plot. But before discussing global linkage analysis, we must first finish dissecting the Scheme 2 partition function.

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

It is useful to subdivide (parse) the Scheme 2 CPFs *Z _{C}* and

*Z*with regard to the resting (R) and activated (A) states of the voltage-sensing J particles, as follows:

_{O}The sub-schemes for the possible configurations of voltage-sensing L and J particles are given by the four (Ca^{2+}-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*, share a single polynomial function, in this case:

_{O}*f*

_{2}(

*λ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*and

_{C}*Z*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

_{O}*g*

_{2}(

*λ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*〉

*. Either parsing route, if connected to a suitable work function, also enables one to measure Δ*

_{V}*W*, as demonstrated shortly using a

_{E}*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^{2+}-binding elements, and is described next.

### The electrical capacitance energy W_{C[q]}

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*of displacing total gating charge per channel (Δ

_{q}*q*) and is equal to the thermodynamic integral:

_{T}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*

_{1}+

*A*

_{3}in Fig. 3 A. The product Δ

*q*

_{T}V_{(+)}is equal to

*A*

_{2}+

*A*

_{3}. The difference between these two areas, −(

*A*

_{2}−

*A*

_{1}), 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*, which was defined in introductory remarks as

_{M}*W*/Δ

_{q}*q*= −

_{T}*W*

_{C}_{[q]}/Δ

*q*. From Eq. 26, we can write the following integral expression for

_{T}*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*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

_{M}*A*

_{6}+

*A*

_{7}−

*A*

_{4}, which represents the value of

*W*/Δ

_{q}*q*and is therefore equal to

_{T}*V*. In turn,

_{M}*V*is also equal to the area

_{M}*A*

_{5}+

*A*

_{6}. Equating the two areas, we obtain

*A*

_{7}=

*A*

_{4}+

*A*

_{5}, which, as anticipated, divides the shaded regions in Fig. 3 B into equal parts. The fact that

*V*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 Δ

_{M}*q*=

_{T}*Q*/

_{max}*N*. Assuming that the value of Δ

*q*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

_{T}*V*(which is independent of

_{M}*N*) can be used to monitor the work function

*W*

_{C}_{[q]}= −Δ

*q*.

_{T}V_{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*, 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 〈

_{M}*q*〉 to

*V*, one obtains an alternative expression for

*V*as a function of the normalized electrical capacitance

_{M}*f*= Δ

_{q}*q*

_{T}^{−1}

*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*, and which in a ligand-binding system is

_{M}*μ*—are intimately related to coefficients in the partition function.

_{M}^{2}Adding the fact that, as we discuss next, there are several practical methods for measuring

*f*that make it easy to compute

_{q}*V*by using Eq. 29, there are compelling reasons to use “mean” in describing

_{M}*V*, although the historical use of “median” is certainly not wrong for the reasons cited. Here, we will refer to

_{M}*V*as the “mean voltage of activation.”

_{M}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*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

_{q}*f*in one sweep, is later demonstrated here through simulation. We note that

_{q}*f*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).

_{q}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*, one can obtain a lower bound (=

_{T}*Δ〈*

^{V}*n*〉

*) for the total number of calcium-binding sites Δ*

_{max}*n*as an alternative to measuring the Hill coefficient.

_{T}### A statistical mechanical formulation of W_{C[q]}

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

^{2+}, and

*Z*

_{V}_{(+)}is analogous except with J and L fixed in their activated positions. The LPFs contain CPFs

*Z*and

_{CR}*Z*, as follows:

_{OA}*Z*

_{V}_{(−)}=

*Z*

_{CR}^{4}and

*Z*

_{V}_{(+)}=

*L*

_{V}_{(+)}(

*J*

_{V}_{(+)}

*DZ*)

_{OA}^{4}. Recalling that

*η*= −

_{L}*kT*ln

*L*and

*η*= −

_{J}*kT*ln

*J*, 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*+ 4Δ

_{L}*q*, the electrical components of

_{J}*η*

_{L}_{,V(+)}and

*η*

_{J}_{,V(+)}exactly cancel Δ

*q*

_{T}V_{(+)}, leaving only the charge-neutral potentials

*η*

_{L}_{,V=0}= Δ

*G*− Δ

_{L}*n*and

_{L}μ*η*

_{J}_{,V=0}= Δ

*G*− Δ

_{J}*n*. In unusual cases where obligatory calcium binding is associated with the activation of L and/or J particles, one would observe sloping asymptotes (

_{J}μ*m*= Δ

*n*+ 4Δ

_{L}*n*) when plotting

_{J}*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*, and recognizing that this is equal to the expectation value 〈

_{CR}*CE*〉

*. The outcome of Eq. 34 solves the earlier problem of not being able to measure*

_{K}*E*solely through conductance Hill analysis, because now the value of 4Δ

*W*can be obtained by subtracting Eq. 34 from Eq. 23a.

_{E}### Limited global analysis

It is reasonable to enquire, given the outcome of Eq. 34, whether a limited version of global analysis can be used to measure Δ*W_{E}* independently of Δ

*W*. The answer is a qualified yes, provided one can maintain the pore in either the closed (

_{C}*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*(∂ln

*Z*

_{L}_{(−)}/∂

*V*)

*and 〈*

_{μ}*q*〉 =

_{O}*kT*(∂ln

*Z*

_{L}_{(+)}/∂

*V*)

*. Using Eq. 12, these can be expanded to:*

_{μ}In BK channels, it is possible to resolve 〈*q _{C}*〉 and 〈

*q*〉 through kinetic means by taking advantage of the fact that pore opening is ∼100 times slower than voltage-sensor activation (Horrigan and Aldrich, 2002).

_{O}Purely thermodynamic methods can also be used to determine 〈*q _{C}*〉 and 〈

*q*〉 from measurements of G-V and Q-V plots, for example, through the relations 〈

_{O}*q*〉 = 〈

_{C}*q*〉 +

*kTd*ln(

*G*−

_{max}*G*)/

*dV*and 〈

*q*〉 = 〈

_{O}*q*〉 +

*kTd*ln

*G*/

*dV*. The expression

*kTd*ln

*G*/

*dV*is known as the mean activation charge displacement 〈

*q*〉, first used in reference to the Shaker K

_{a}^{+}channel (Sigg and Bezanilla, 1997), where it appears to follow the course of an “upside-down” Q-V curve. In BK, measurements of 〈

*q*〉 (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

_{a}_{o}), compatible with that of Δ

*q*. In contrast, the limiting slope observed in Shaker is substantially larger (12–13 e

_{L}_{o}; see Islas and Sigworth, 1999)—equal in value to Δ

*q*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.

_{T}A different but related method uses the slope of the conductance Hill to determine 〈*q _{C}*〉 and 〈

*q*〉. Starting from the identity

_{O}*W*

_{H}_{[g]}=

*kT*ln(

*Z*

_{L}_{(+)}/

*Z*

_{L}_{(−)}), we derive

*dW*

_{H}_{[g]}/

*dV*= 〈

*q*〉 − 〈

_{O}*q*〉 (see also Conti, 1986), and recognizing that 〈

_{C}*q*〉 = (1 −

*P*)〈

_{O}*q*〉 +

_{C}*P*〈

_{O}*q*〉, we obtain 〈

_{O}*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*for the simple reason that the J–K interaction does not directly involve the pore.

_{E}### Model-independent application of linkage analysis

Linkage plots yield model-independent information about channel gating, starting with whether the channel is part of an “accessible” allosteric network. We have just discussed how 〈*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*) and positive (

_{T}*m*

_{V}_{(+)}= Δ

*q*) asymptotes.

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

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

Specifically, the voltage sensor and pore are both voltage sensitive (Δ*q _{L}*, Δ

*q*≠ 0), whereas the calcium-binding domain demonstrates no voltage sensitivity (Δ

_{J}*q*= 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.”

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

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

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

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

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

This last point is the weakest in terms of experimental evidence supporting it. There is no a priori reason why allosteric interactions should be independent of environmental influence. It has been postulated, for example, that interactions between voltage-sensing particles in the voltage sensor might also be voltage sensitive (Pantazis et al., 2010), and interactions between the gating ring and the voltage sensor appear to be mediated by internal Mg^{2+} (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**,

**) and internal allosteric factors (**

*K***,**

*F***), as indicated by bold type. The linkage diagram of Scheme 3 in the parsing order**

*G**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

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

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

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

The test model (Scheme 4) used in kinetic simulations is shown in Fig. 5. The scheme consists of a central pore (L) surrounded by four identical regulatory subunits, with each subunit containing two voltage-sensing particles (J1 and J2) and two calcium-sensitive particles (K_{1} and K_{2}), for a total of 17 interacting gating particles. The particles L, J_{2}, K_{1}, and K_{2} form an allosteric loop involving the core interactions *C*, *D*, and *E*, and an internal interaction *F* coupling K_{1} and K_{2}. The J_{1} particle does not participate in the loop but influences its neighboring J_{2} 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_{2}→(J_{1}, K_{2})→K_{1} is given by:

from which we derive, as usual, the first-order equation *Z* = *Z _{C}*

^{4}+

*LZ*

_{O}^{4}. Expressions for the CPFs

*Z*and

_{C}*Z*are as follows:

_{O}The second-order CPFs *Z _{CR}*,

*Z*,

_{CA}*Z*, and

_{OR}*Z*are functions of

_{OA}*K*

_{1}and

*K*

_{2}, given by

*Z*=

_{CR}*f*

_{2}(

*K*

_{1},

*K*

_{2}) = (1 +

*K*

_{1}) +

*K*

_{2}(1 +

*K*

_{1}

*F*),

*Z*=

_{CA}*f*

_{2}(

*K*

_{1}

*E*,

*K*

_{2}),

*Z*=

_{OR}*f*

_{2}(

*K*

_{1},

*K*

_{2}

*C*), and

*Z*=

_{OA}*f*

_{2}(

*K*

_{1}

*E*,

*K*

_{2}

*C*). The mean single-channel gating charge displacement 〈

*q*〉 ranges from 0 to Δ

*q*= Δ

_{T}*q*+ 4(Δ

_{L}*q*

_{J}_{1}+ Δ

*q*

_{J}_{2}) and is calculated from 〈

*q*〉 = 〈

*l*〉Δ

*q*+ 〈

_{L}*j*

_{1}〉Δ

*q*

_{J}_{1}+ 〈

*j*

_{2}〉Δ

*q*

_{J}_{2}, where the voltage-dependent particle activation curves are derived as follows:

The average single-channel conductance is 〈*g*〉 = 〈*l*〉*g _{L}*.

### Linkage analysis of simulated data using Scheme 4

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

^{2+}]

_{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*, and three out of five averaged outcomes found to be statistically indistinguishable from their true values.

_{E}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*= 2.62 e

_{T}_{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*, and

_{eq}*g*. From the definition

_{leak}*W*

_{H}_{[g]}=

*kT*ln

*G*−

*kT*ln(

*G*−

_{max}*G*), we recognize that the

*kT*ln

*G*defines the

*V*

_{(−)}asymptote of

*W*

_{H}_{[g]}and is logarithmically sensitive to the floor, whereas

*kT*ln(

*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*and variance

_{L}P_{O}*Var*=

_{G}*Ng*(1 −

_{L}P_{O}*P*) through the relation S/N =

_{O}*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

*kT*ln

*G*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*< 1/

_{O}*N*. Therefore, with a cumulative

*N*of 10

^{5}(= 10

^{2}traces × 10

^{3}channels), noise effects became evident for

*P*≤ 10

_{O}^{−5}, as seen in Fig. 9. It is remarkable that useful data could be derived for

*P*smaller than the predicted limit, down to

_{O}*P*= 10

_{O}^{−8}. 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

^{5}, noise is increased when 1 −

*P*= 1/

_{O}*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*, it was least affected by floor and ceiling effects and was used to fit the slope Δ

_{O}*q*. The average fitted value of Δ

_{L}*q*differed only slightly (+1.6%) from the given value of 0.3 e

_{L}_{o}. To determine interaction energies Δ

*W*and Δ

_{C}*W*, the relative heights of the other three “noisy” asymptotes were estimated through a process of alternatively fitting and manually adjusting the patch-specific variables (

_{D}*N*,

*V*, and

_{eq}*g*) so that they visually paralleled the (

_{leak}*V*

_{(−)},

*μ*

_{(+)}) asymptote. The value of Δ

*W*(averaged from measurements taken at

_{C}*V*

_{(−)}and

*V*

_{(+)}) was 1.3% less than the given value of −53 meV. Measuring Δ

*W*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

_{D}^{2+}plot was excluded from the analysis.

In *W _{C}*

_{[q]}experiments (Fig. 10), the average value of the outcome parameter Δ

*W*+ Δ

_{C}*W*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

_{E}*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

^{2+}]

_{i}. It was quantified by simulating the model in the absence of Nyquist noise, and the

*Δ*

^{V}*V*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 Δ

_{M}*q*measurements, accounted for most of the error in Δ

_{T}*W*+ Δ

_{C}*W*.

_{E}### Hill analysis performed on markers of voltage- and Ca^{2+}-sensor activation

Hill analysis is not limited to G-V plots, provided that markers of activation can be obtained for particles other than the pore. Energy Hill plots for each of the four subparticles (J_{1}, J_{2}, K_{1}, and K_{2}) in Scheme 4 are shown in Fig. 11. The work functions for the Ca^{2+}-sensitive K particles are plotted on *μ* and log[Ca^{2+}]_{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_{1} providing the best resolution because of its small intrinsic charge and relatively large interaction energy with J_{2}. From an experimental standpoint, the Hill plot for J_{2} 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_{1} and K_{2} (Δ*W _{F}* > 0).

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

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

As was done previously for K particles in Scheme 1, multiplicative constants *B*_{1}, *B*_{2}, *B*_{3}, and *B*_{4}, 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*/4 was additionally offset by the amount −Δ

_{L}*W*= −(Δ

_{B}*W*

_{B}_{1}+ Δ

*W*

_{B}_{2}+ Δ

*W*

_{B}_{3}+ Δ

*W*

_{B}_{4})/4 to maintain the original free energy change Δ

*G*between closed and open states. Under these conditions, of the six configurational pore states (which include the cis- and trans-configurations of the

_{L}*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*and −2Δ

_{B}*W*; Fig. 12 B) and brief occupancy times. These attributes of the interacting tetrameric pore are found in the Scheme 5 partition function:

_{B}where *Z _{C}* and

*Z*remain unchanged from Scheme 4 (Eq. 38). It is readily apparent that for large values of

_{O}*B*= exp(−Δ

*W*/

_{B}*kT*), the partition function reverts back to Scheme 4 with its binary pore equilibrium constant

*L*=

*P*

^{4}. 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*/4, 〈

_{L}*p*〉 = ∂ln

*Z*/∂ln

*P*, and

*η*=

_{P}*η*/4. The results of global capacitance (

_{L}*W*

_{C}_{[q]}) analysis were practically unchanged from Scheme 4 because interactions between voltage- and Ca

^{2+}-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*+ 2Δ

_{D}*W*). The value of −2Δ

_{B}*W*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 Δ

_{B}*W*in native BK channels would be considerable, considering the large potentials needed to reach “true” asymptotic behavior and the need to resolve subconductance levels.

_{B}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*. The loss of uniform incremental values in these parameters influences the measurement of self-interaction energy, as demonstrated by the new linkage relation:

_{p}g_{L}The second term on the right of Eq. 42 is an artifact of the nonuniform spacing of subconductance levels and vanishes if *x*_{1} = 0.25 and *x*_{3} = 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: *kT*ln(1/B) and *kT*ln(*B*^{3}/*B*_{1}*B*_{2}), the difference of which is *kT*ln(*B*_{3}*B*_{4}) = −(Δ*W _{B}*

_{3}+ Δ

*W*

_{B}_{4}). These expressions stem from the constraints imposed on self-interaction energies to achieve a “balanced” equilibrium between closed and open states.

### Dual-modality particles

Constraint 1 for BK gating required all gating particles to be unimodal, i.e., each responsive to only one external force. However, it is not hard to imagine circumstances in which this would be false. For example, a Ca^{2+}-binding site might be displaced from the internal solution by a fraction *d* of the membrane potential, leading to a gating charge displacement of 2*d* e_{o}. Another mechanism for dual sensitivity is the obligatory binding of Ca^{2+} during pore opening. As unlikely as this seems, Piskorowski and Aldrich (2002) have reported Ca^{2+} 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^{2+} upon opening was performed by specifying the particle function *μ _{L}* = Δ

*G*+ −Δ

_{L}*q*− Δ

_{L}V*n*, where Δ

_{L}μ*n*was chosen to be 0.4.

_{L}^{3}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).

### Modality-sensitive allosteric factors

Constraint 4 implied that allosteric energies are insensitive to external forces. There is reason to believe that this might not be the case in BK. For example, Mg^{2+} 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^{2+} 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*. A voltage-sensitive

_{G}V*G*= exp(−Δ

*G*/

_{G}*kT*) would contribute an amount (∂ln

*Z*/∂ln

*G*)Δ

*q*of gating charge to the Q-V curve. The effects on the Q-V and Hill plots from a small charge (Δ

_{G}*q*= 0.2 e

_{G}_{o}) are shown in Fig. 14 (A and B). Because interactions between voltage- and Ca

^{2+}-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*grew larger by the amount 4Δ

_{T}*q*. On the other hand, the

_{G}*V*asymptotes in the local Hill plot for J

_{2}diverged because of the addition of Δ

*q*to the slope of the upper asymptote (Fig. 14 B). The interaction energy −(Δ

_{G}*W*+ Δ

_{G}*W*) could still theoretically be obtained by evaluating the height differences between asymptotes evaluated at

_{D}*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

_{2}and K

_{1}). Here, a more dramatic change to the Q-

*V*curve occurred (Fig. 13 C). At very low [Ca

^{2+}]

_{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*being shifted immeasurably far to the right on the

_{M}*V*axis. With increasing [Ca

^{2+}]

_{i}, The

*E*charge component (∂ln

*Z*/∂ln

*E*)Δ

*q*entered the Q-V plot from the right side, with

_{E}*V*eventually reaching a well-defined limiting value specified by

_{M}*V*

_{Mμ}_{(+)}. The global energy of interaction −4(Δ

*W*+ Δ

_{C}*W*) can be deduced from the area

_{E}*A*

_{2}−

*A*

_{1}shown in Fig. 14 C. From an experimental standpoint, obtaining the correct result would require that one recognize that Δ

*q*, the total charge movement per channel, has become a function of [Ca

_{T}^{2+}]

_{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.

### Intersubunit interactions

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

## DISCUSSION

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*, 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

_{H}*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.

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

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

### Sources of error in electrical capacitance energy measurements

The major source of error in *W _{C}*

_{[q]}measured in voltage-ramp simulations was systematic, caused by a Ca

^{2+}-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*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

_{M}*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

^{6}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*, 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 (=

_{max}*mC*) 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

_{patch}^{−4}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*(assuming Δ

_{J}*C*is voltage independent) would obtain a second-order voltage term:

_{J}*η*= Δ

_{J}*G*− Δ

_{J}*q*− Δ

_{oJ}V*C*(

_{J}*V*−

*V*)

_{o}*V*. This could theoretically have an effect on baselines and limits of integration when analyzing ramp-generated

*C*to measure

_{q}*V*. 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).

_{M}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*= 2.62 e

_{T}_{o}has been determined for BK (Horrigan and Aldrich, 2002), an independent determination of Δ

*q*using charge-per-channel methods has yet to be published. Such experiments (see Fig. 8) require that one measure gating currents and

_{T}*N*in the same preparation (Aggarwal and MacKinnon, 1996; Seoh et al., 1996).

### Sources of error in conductance Hill energy measurements

The critical step in *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*< 1/

_{O}*N*for the floor asymptote, and 1 −

*P*= 1/

_{O}*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^{−7} (see Fig. 9), requiring at the very least a 24-bit A/D converter (10^{7} ≈ 2^{24}). 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*and baseline

_{eq}*I*of the ionic current, and then applies a logarithmic converter before digitizing. Using a repeating ramp protocol, manual or semi-automated adjustment of inputs

_{base}*G*,

_{max}*V*, and

_{eq}*I*could be performed to achieve an optimal shape of the Hill curve in real time.

_{base}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*. 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).

_{O}### Relation to the HA model

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

From measurements of the Ca^{2+} dependence of *P _{O}* at very negative potentials, the authors computed a variable they called

*R*

_{O}that is equivalent to

*Z*/

_{OR}*Z*, and by plotting log

_{CR}*R*

_{O}versus [Ca

^{2+}] obtained log〈

**〉. They also obtained the intrinsic pore charge Δ**

*C**q*by fitting the negative asymptote of log(

_{L}*P*). They recognized the importance of using log(

_{O}*P*), 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

_{O}*P*(using an all-points histogram of single-channel currents to measure

_{O}*NP*and dividing by

_{O}*N*≈

*G*/

_{max}*g*) did not lend itself to accurate measurements at the ceiling of the G-V curve (a requirement for defining the positive

_{L}*V*asymptote of the conductance Hill plot). They cited this as a difficulty in estimating the value of

**, defined in linkage terms as the height difference in the two asymptotes of**

*D**W*

_{H}_{[g]}versus

*V*. Instead, an alternative method for estimating

**using the product of Δ**

*D**q*and the midpoint voltages of 〈

_{J}*q*〉 and 〈

_{C}*q*〉 was proposed specifically for the HA model (Scheme 2). The proof for the more general Scheme 3 goes as follows: recognizing that 〈

_{O}*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*〉 into its baseline (Δ

_{O}*q*) and variable (〈

_{L}*q*〉 − Δ

_{O}*q*) components, we subtract

_{L}*Δ*

^{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*〉 and 〈

_{C}*q*〉 − Δ

_{O}*q*, respectively. Recalling that

_{L}*Δ(*

^{V}*W*

_{H}_{[g]}+

*η*) = −4

_{L}*kT*ln(

**), and defining**

*D**Δ*

^{L}*V*=

_{M}*V*

_{M}_{(O)}−

*V*

_{M}_{(C)}, we obtain the desired result: Δ

*q*

_{J }*Δ*

^{L}*V*= 4Δ

_{M}*W*(see Fig. 15). Assuming one can measure both 〈

_{D}*q*〉 and 〈

_{C}*q*〉—for example, by integrating the fast component of gating charge after prepulsing to either

_{O}*V*

_{(−)}or

*V*

_{(+)}—then Δ

*W*could be calculated without resorting to logarithmic (Hill) analysis of the G-V curve. The value of Δ

_{D}*q*, which is necessary to complete the calculation, can theoretically be obtained by subtracting Δ

_{J}*q*from Δ

_{L}*q*, both of which can be independently measured. The authors acquired 〈

_{T}*q*〉 indirectly through the relation 〈

_{O}*q*〉 = 〈

_{O}*q*〉 + 〈

_{a}*q*〉, where 〈

*q*〉 =

_{a}*kTd*ln

*G*/

*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(

**) compared with the fitted value from log(**

*D**P*)-

_{O}*V*curves.

The allosteric factor ** E** was estimated in similar fashion by measuring shifts in 〈

*q*〉 and 〈

_{C}*q*〉 in response to saturating changes in [Ca

_{O}^{2+}]

_{i}. Because the influence of pore activation has been removed from both 〈

*q*〉 and 〈

_{C}*q*〉, 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 Δ

_{O}*W*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

_{G}*P*was very low (∼10

_{O}^{−7}), is reassuring for the eventual success of linkage methods described here. It will be of interest to see if, despite its potential limitations, the slow voltage ramp can overcome the inability of the

*NP*histogram to generate positive voltage (“ceiling”) asymptotes in the conductance Hill plot.

_{O}### Deconstructing other thermodynamic methods of analysis

Methods for dealing with equilibrium plots are prevalent in the ion channel literature. The tools used most often are the Hill and Boltzmann equations, which can provide very good fits to experimental data. Confusion arises from different uses of the “Hill” label. The traditional Hill plot is defined as ln(*M*/(*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 ^{n}*, 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

^{n}*K*/(1 +

^{n}*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}*n*is the fitted parameter, named

*n*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 Δ

_{H}*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*and the expression

_{M}*W*

_{C}_{[q]}= −Δ

*q*, as described here and in Chowdhury and Chanda (2012).

_{T}V_{M}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*is often used as a stand-in for the more correct

_{1/2}*V*. In a somewhat confusing twist of nomenclature, the

_{M}*V*(the voltage for which the normalized Q-V is 0.5) is equal to the median value of the capacitance distribution function

_{1/2}*f*(see accompanying article by Chowdhury and Chanda in this issue), whereas

_{q}*V*, 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,

_{M}*V*has no special thermodynamic significance (in the sense of describing an energy or displacement), but it may coincide with

_{1/2}*V*for symmetrical Q-V plots.

_{M}Similarly, the product Δ*q _{B}V_{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*〉, where 〈

_{a}*q*〉 =

_{a}*kTd*ln

*G*/

*dV*is the mean activation charge displacement. Therefore, gating current measurements are not required to measure

*W*

_{C}_{[q]}, because even Δ

*q*can be obtained from the limiting value of activation charge displacement: Δ

_{T}*q*= 〈

_{T}*q*〉

_{a}

_{V}_{(−)}. However, the logarithmic dependence of 〈

*q*〉 on

_{a}*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

^{2}ln

*G*/

*dV*

^{2})

*dV*. This contrasts with the traditional Boltzmann method of analysis, where Δ

*q*

_{B}V_{1/2}obtained by Boltzmann fitting is numerically equal to −

*kT*ln(

*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*

_{L}V_{M}_{[g]}, where

*V*

_{M}_{[g]}= ∫

*Vf*and

_{g}dV*f*=

_{g}*dP*/d

_{O}*V*. The variable

*V*

_{M}_{[g]}could be referred to as the “mean” activation voltage of conduction, and

*f*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 Δ

_{g}*q*Δ

_{L}^{μ}*V*

_{M}_{[g]}, as a kind of “global analysis” of the local K–L interaction, is equal to −4Δ

*W*. This attempt at deriving the strength of K–L linkage (or any perturbation affecting pore opening, say from a mutation) by measuring shifts in

_{C}*V*

_{M}_{[g]}is flawed, because, unlike

*Q*=

*NkT*(∂ln

*Z*/∂ln

*V*), the conductance

*G*= (∂ln

*Z*/∂ln

*L*)

*G*cannot be expressed as the voltage derivative of

_{max}*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*− Δ

_{O}dV*q*

_{L}V_{(+)}, 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*(∂ln

*Z*/∂

*η*) points to the underlying problem, which is that

_{L}*η*cannot be varied independently of

_{L}*η*because both particle potentials are functions of the control parameter

_{J}*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*(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.

_{E}To summarize this section, we have excluded certain candidates from a list of valid work functions. These include: (a) the product Δ*q _{B}V_{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.

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

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

*H*− Δ

_{K}*S*, where Δ

_{K}T*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*. The steepness of the curve is determined by Δ

_{K}*S*. Rearranging

_{K}*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

^{55}, which is equivalent to an entropy change of Δ

*S*= 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 Δ

_{K}*H*= Δ

_{K}*S*≈ −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.

_{K}T_{o}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*can be measured from the rise of the sigmoidal component of the conductance Hill energy

_{C}*W*

_{H}_{[g]}plotted versus

*T*. If there is also an entropic component Δ

*S*, it would show as a difference in the slopes of the positive and negative asymptotes equal to Δ

_{C}*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*evaluated at

_{Ko}*T*(Privalov and Khechinashvili, 1974), which would generate a measurable nonlinearity in the particle potential

_{o}*η*of a T-sensor. The value of Δ

_{K}*C*appears to be fairly constant (Baldwin, 1986), leading to

_{p}*η*= (Δ

_{K}*C*− Δ

_{p}*S*)(

_{Ko}*T*−

*T*) − Δ

_{o}*C*ln(

_{p}T*T*/

*T*), which is derived by integrating the thermodynamic identities

_{o}*C*= (∂

_{p}*H*/∂

*T*)

*=*

_{P}*T*(∂

*S*/∂

*T*)

*around*

_{P}*T*= Δ

_{o}*H*/Δ

_{Ko}*S*, evaluated for both native and denatured states (Clapham and Miller, 2011).

_{Ko}As a consequence of the change in heat capacitance with melting, dramatic changes are seen in *W _{H}*

_{[g]}related to the possibility that

*η*crosses zero twice, theoretically generating two melting temperatures separated approximately by 2Δ

_{K}*T*= 2Δ

*H*/Δ

_{Ko}*C*, where Δ

_{p}*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*, the sensor would be hot-activated (negative slope for

_{o}*η*), 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.

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

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

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

On the other hand, there are advantages to considering thermodynamic models and the linkage techniques used to analyze them. The measured quantities in linkage analysis are energy and displacement, which are used to characterize gating networks consisting of regulatory particles. Kinetic analysis, on the other hand, measures rate constants and displacements for a particular scheme of connected microscopic states. The partition function, which defines the thermodynamic properties of a channel, can often be expressed as a simple polynomial, even for complex allosteric models such as Schemes 4 and 5 that possess >10^{3} 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*〉 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.

_{a}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.

## APPENDIX

### Linkage diagrams

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)\u2192L(KC)$⇒(1 + *K*) + *L*(1 + *KC*). If *C* = 1, implying independence between L and K, then the CPF becomes (1 + *K*)(1 + *L*), which can also be written in linkage notation as (K)(L). Using these rules and the linkage diagrams found in the text, the partition function for Scheme 1 is seen to be *Z* = (1 + *K*)^{4} + *L*(1 + *KC*)^{4}, and for Scheme 2, it is *Z* = [(1 + *K*) + *J*(1 + *KE*)]^{4} + *L*[(1 + *KC*) + *JD*(1 + *KCE*)]^{4}. Linkage diagrams are not unique. For example, $(K)\u2192L\u2009(KC)$ is equivalent to $(L)\u2192K\u2009(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^{2+} sensors *K*_{1}, *K*_{2}, *K*_{3}… can be represented by *K*_{[F]}, where ** F** = {

*F*

_{1},

*F*

_{2},

*F*

_{3}…} 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*

_{1}

*K*

_{2}

*K*

_{3}…

*F*

_{1}

*F*

_{2}

*F*

_{3}…). 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.

## Acknowledgments

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.

## References

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

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

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