## Abstract

Ion channels are membrane-bound enzymes whose catalytic sites are ion-conducting pores that open and close (gate) in response to specific environmental stimuli. Ion channels are important contributors to cell signaling and homeostasis. Our current understanding of gating is the product of 60 plus years of voltage-clamp recording augmented by intervention in the form of environmental, chemical, and mutational perturbations. The need for good phenomenological models of gating has evolved in parallel with the sophistication of experimental technique. The goal of modeling is to develop realistic schemes that not only describe data, but also accurately reflect mechanisms of action. This review covers three areas that have contributed to the understanding of ion channels: traditional Eyring kinetic theory, molecular dynamics analysis, and statistical thermodynamics. Although the primary emphasis is on voltage-dependent channels, the methods discussed here are easily generalized to other stimuli and could be applied to any ion channel and indeed any macromolecule.

Ion channels are membrane-bound enzymes whose catalytic sites are ion-conducting pores that open and close (gate) in response to specific environmental stimuli (voltage, ligand concentration, membrane tension, temperature, etc.; Hille, 2001). These stimuli activate highly specialized regulatory domains (an example is the voltage-sensing domain [VSD] found in many cation-selective channels) that are coupled to the pore gate through as yet incompletely understood mechanisms. Ion channels are important contributors to cell signaling and homeostasis and are strictly necessary for electrical conduction in nerve and muscle tissue. Our current understanding of channel gating is the product of over 60 years of voltage-clamp recording augmented by experimental intervention in the form of environmental, chemical, and mutational perturbations. Macroscopic ionic or capacitive gating currents reflecting the collective behavior of many channels can be interpreted with the help of descriptive (phenomenological) models that range from the standard two-state Boltzmann fit to more sophisticated analysis using multistate kinetic schemes (Hodgkin and Huxley, 1952; Vandenberg and Bezanilla, 1991; Zagotta et al., 1994; Schoppa and Sigworth, 1998b; Horrigan and Aldrich, 2002). Elementary gating events are resolvable using fluctuation analysis (Neher and Stevens, 1977; Sigworth, 1980, 1981; Conti and Stühmer, 1989; Sigg et al., 1994) or by measuring single-channel ionic currents (Colquhoun and Hawkes, 1995a). More recently, fluorometry and other optical techniques have been added to the repertoire of electrophysiological recordings (Mannuzzu et al., 1996; Cha and Bezanilla, 1997; Perozo et al., 1999). The tradition of modeling electrophysiological data dates back to the early 50s, when Hodgkin and Huxley published their celebrated 1952 papers culminating with a general description of the squid giant axon action potential (Hodgkin and Huxley, 1952). Alan Hodgkin writes about this time period (Hodgkin, 1976): “As soon as we began to think about molecular mechanisms it became clear that the molecular data would by itself yield only very general information about the class of system likely to be involved. So we settled for the more pedestrian aim of finding a simple set of mathematical equations which might plausibly represent the movement of electrically charged gating particles.”

What began as a quest to understand the physics of axon excitability (Hodgkin and Huxley had originally favored a carrier scheme) evolved into the slightly less ambitious objective of computing the time course of the action potential using mathematically defined units of activation, each containing four independent gating particles, thus anticipating the four-subunit structural motif now known to describe most ion channels. A plan to improve the fitting of voltage-clamp records by adding more gating particles was apparently tempered by the laborious necessity of performing numerical integration with a hand-cranked calculator: “Better agreement might have been obtained with a fifth or sixth power, but the improvement was not considered to be worth the additional complication” (Hodgkin and Huxley, 1952).

The past 60 years have seen tremendous advances in technique, including the ability to clone, mutate, and insert channels of any species into heterologous expression systems and the determination of x-ray crystal structures (Doyle et al., 1998; Long et al., 2005) in selected channels (and faster computers!). This has put us on the path to a residue-level understanding of ion channel structure and function. The need for good phenomenological models has evolved commensurately. More than two centuries after Luigi Galvani (1737–1798) formulated his “electrical hypothesis” of muscle contraction, we are approaching the threshold of a molecular and mechanistic description of membrane excitability. The goal is to develop gating schemes that do more than describe data: they also accurately reflect structures and mechanisms of action. This is an opportune time, then, to reflect on various kinetic and equilibrium models that have been used in the past, note recent advances, and comment on possible directions for the future of ion channel modeling.

The aim of this review is to explore common ground in the phenomenology of three fields that have contributed to the study of ion channels: traditional Eyring kinetic theory, statistical thermodynamics, and, most recently, molecular dynamics (MD) analysis. Much of the text will be devoted to the application of these fields as expressed by the following quantities: the transition rate constant (α), the potential of mean force (*W* [PMF]), and the partition function (*Z*). The emphasis will be on activated (barrier) processes with slow (millisecond) relaxation rates (readers wishing instead to review low-barrier electrodiffusion may start with these references: Bowman and Baglioni, 1984; Im and Roux, 2002). An essential component of the stated goal is a shared terminology. In view of the central role played by the voltage clamp in ion channel experimentation, the discussion of kinetics will use electrical analogues to the usual mechanical variables; for example, charge (*q*) ↔ distance, voltage (*V*) ↔ force, inductance (*L*) ↔ mass, and so on. However, recognizing that ion channels demonstrate tremendous variety in their responses to external stimuli, a subsequent section on thermodynamics is developed around canonical displacement-force pairs (ξ, Ω), thereby broadening the discussion from voltage-dependent (*q*, *V*) channels to ligand-gated (*n*, μ) and other types of channels. The scope of ion channel physiology is vast, and this review cannot hope to cover every facet of its phenomenology. As a former experimentalist turned itinerant modeler, the author does not claim a practicing knowledge in many areas touched upon here; in these an experimental viewpoint is adopted. Hopefully, the reader will be able to fill in any deficiencies from references spanning 80 years of theoretical and experimental work.

### Ion channel kinetics

Kinetic modeling of ion channels began with the Hodgkin–Huxley (HH) scheme (Hodgkin and Huxley, 1952) and continues unabated today with the widespread use of multistate Markov models and single-channel analysis. The “pedestrian aim” expressed by Hodgkin and Huxley, namely to simulate the shape and threshold nature of stimulated action potentials through simple kinetic equations, succeeded brilliantly, setting the stage for countless papers applying the HH formulism to model excitability in a variety of tissues. The building block of the HH model is the gating particle, which transitions between resting and activated states in a voltage-dependent manner. Because ion channels are composed of modular domains, each with its own specialization and distinct evolutionary origin (Schreiber et al., 1999; Anderson and Greenberg, 2001; Murata et al., 2005; Sasaki et al., 2006), the gating particle concept remains a viable one, although gating schemes have grown in complexity. The opening of a potassium (K^{+}) channel in the original HH scheme required four randomly fluctuating “n” particles to be simultaneously active (the sodium channel description was similar but used three activating “m” particles and one inactivating “h” particle). K^{+} conduction in the HH model is therefore proportional to *n*^{4}, with particle kinetics satisfying *n*_{∞} = α/λ are functions of the voltage-dependent unidirectional rate constants α and β. The HH phenomenological equations are consistent with the notion that gating particles undergo transitions between macrostates (Fig. 1 A). The invention of the patch clamp (Neher and Sakmann, 1976) allowed single-channel experiments to confirm the existence of stochastic transitions between discrete conductance levels as predicted by the HH model.

The physics of particle transitions concerns the study of chemical escape from a metastable state, known as reaction rate theory, which has a long and interesting history (reviewed by Hänggi and Borkovec, 1990). The particular flavor of reaction rate theory best known to the ion channel community is the transition state theory (TST) of Eyring, which is premised on the somewhat tenuous assertion that reactants (R) rapidly thermalize with their surroundings until they reach the separatrix of the transition barrier (b), whereupon they inexorably turn to product (P). A one-dimensional bistable landscape serving as a model of unimolecular gating particle kinetics is shown in Fig. 1 A. The TST rate constant for the forward reaction is*kT*/*h* equals the reactive frequency of barrier crossing, the activation energy is Δ*E** = *E _{b}* −

*E*, and the

_{R}*Z*s are statistical weights (partition functions) for the translational, vibrational, and rotational degrees of freedom of the reactant and barrier states (Eyring, 1935a,b). Pressure-volume work added to Δ

*E** yields the activation enthalpy Δ

*H** = Δ

*E** + Δ

*v**

*P*. The barrier partition function

*Z*′

*is primed to indicate that it is missing the unstable forward translational component*

_{b}*Z*= (2π

_{tran}*mkT*)

^{1/2}/

*h*per unit length that, when multiplied by the mean forward velocity (

*kT*/2π

*m*)

^{1/2}, is the source of the

*kT*/

*h*prefactor (Eyring, 1935b). The stable degree of freedom in

*Z*matching

_{R}*Z*is the vibrational function

_{tran}*Z*= 2π

_{vib}*kT*/

*h*ω

*, where ω*

_{R}*is the oscillator frequency of the reactant state. Separating*

_{R}*Z*from

_{vib}*Z*reestablishes term balance in the ratio of statistical weights, yielding the more meaningful prefactor (ω

_{R}*/2π)(*

_{R}*Z*′

*/*

_{b}*Z*′

*). The rate constant is further simplified by moving the*

_{R}*Z*ratio into the exponent, resulting in

*G** = Δ

*H** − Δ

*S**

*T*is the Gibbs activation energy, whose entropy Δ

*S** =

*k*ln(

*Z*′

*/*

_{b}*Z*′

*) contains the nonreactive degrees of freedom. Contributions from external forces such as voltage*

_{R}*V*or chemical potential μ add new terms to Δ

*G** and generate the expanded free energy change Δ

*W** = Δ

*G** − Δ

*q**

*V*− Δ

*n**μ −…, where Δ

*q** (charge) and Δ

*n** (binding number) are the canonical activation displacements to

*V*and μ and other force-displacement pairs may be added.

The major failing of the TST (apart from Eyring’s misleading emphasis on the universal prefactor *kT*/*h*, a mere accounting trick) is that it does not properly account for the influence of the solvent or heat bath. Solvent interactions contribute not only additional degrees of freedom to Δ*S**, but also generate fluctuations in the reacting coordinate that might reduce the rate of reaction (because not all barrier crossings lead to product). Eyring recognized the latter as a deficiency in his theory and suggested an empirical fix by multiplying *k _{TST}* by a “transmission coefficient” κ. With these considerations, Eyring’s rate constant for the one-dimensional barrier assumes its final form:

Despite its theoretical shortcomings, Eyring’s theory has been the mainstay of rate constant phenomenology in the ion channel literature. Voltage-dependent rate constants are usually expressed as variants of the following Eyring-like equations:*q* = Δ*q*_{α} + Δ*q*_{β}. The preexponential factors in Eq. 2 (a and b) are related to Eyring’s theory through α_{o} = κ(ω* _{R}*/2π)exp(−Δ

*G**/

*kT*) and β

_{o}= κ(ω

*/2π)exp[−(Δ*

_{P}*G** − Δ

*G*)/

*kT*], where Δ

*G*=

*G*−

_{P}*G*.

_{R}Channel gating is typically modeled by a network of discrete states, with each state connected to its neighbors through unimolecular rate constants. Such a network may be referred to as a discrete-state Markov (DSM) model. The word “Markovian” describes the absence of “memory,” in which predictions of future events depend solely on the present state. The number of states in DSM gating schemes is typically small (the HH model of the K^{+} channel consists of a linear arrangement of *n* = 5 states, and many empirical schemes do not exceed 10 states, although with allosteric models *n* may rise considerably). The existence of discrete closed-state transitions is supported by fluctuation analysis of gating currents (Conti and Stühmer, 1989; Sigg et al., 1994). Alternative gating models (Millhauser et al., 1988; Liebovitch, 1989) characterized by a “rough” energy landscape designed to simulate power-law kinetics seen in low-temperature protein motions (Frauenfelder et al., 1991) have not gained traction, primarily because the multiexponential kinetics predicted by DSM models adequately describe most electrophysiological data (McManus et al., 1988). The kinetics of an *n*-state DSM model are obtained by solving the master equation d**p**/d*t* = **pA**, where **p**(*t*) is the 1 × *n* vector of state probabilities and **A** is an *n* × *n* array of rate constants α* _{ij}* known as the “Q-matrix” (Colquhoun and Hawkes, 1995b). Knowing

**A**, one can derive the nonstationary mean value and flux for any observable

*a*through the vector relations 〈

*a*〉 = 〈

**p**|

**a**〉 and

**pA**|

**a**〉, where

**a**is the

*n*× 1 state vector of

*a*values, and

**p**(

*t*) is solved as an initial value problem through numerical integration or eigenvector decomposition. The latter method yields

**u**

_{r}and

**v**

_{r}are the left and right eigenvectors of

**A**and the corresponding eigenvalues −λ

_{r}are the macroscopic decay rates. In a closed (gating) system satisfying detailed balance, all λ

_{r}are distinct, real, and positive, except for the zeroth eigenvalue λ

_{0}= 0, whose left eigenvector

**u**

_{0}is the equilibrium probability distribution

**p**(∞). Eigenvalues for the simple HH K

^{+}channel scheme are integer multiples of λ = α + β. In contrast, an open system (pore, transporter) may decay with damped oscillations (Re[λ

_{r}] ≥ 0) toward a nonequilibrium steady-state (Schnakenberg, 1976). An important element of Q-matrix theory is kinetic analysis of single-channel data, but this subject is too vast to cover here (reviewed by Qin, 2007).

Although DSM models have been immensely helpful in studying gating at the network level, the physical origin of the rate constant, specifically the basis for the Eyring variables κ, Δ*q*, ω, and Δ*G**, remains poorly understood. One obstacle to experimentally determining rate constants is the difficulty in measuring single transition rates without the aid of empirical models. An important step forward was achieved when crystal structures of ion channels were made available (Doyle et al., 1998; Long et al., 2005), providing the starting point for all-atom MD simulations and the theoretical ability to compute rate constants in silico. Although most of the MD work on ion channels has focused on ion permeation in pores (Roux and Karplus, 1991; Allen et al., 1999; Shrivastava and Sansom, 2000; Bernèche and Roux, 2001), long-time scale trajectories of voltage sensor activation have recently been published (Delemotte et al., 2011; Jensen et al., 2012), setting the stage for detailed MD studies of gating.

Because of the large time scale difference between the step size (approximately femtoseconds) of MD integration and intervals between ion hops in pores (nanoseconds to microseconds) or the relaxation time of a gating particle (approximately milliseconds), special techniques are required to compute rate constants from MD trajectories. The progenitor of these techniques is the reactive flux method of Chandler (1978), which utilizes as the source of calculation an ensemble of short MD trajectories originating at the barrier peak. With the right choice of reaction coordinate *q* and provided the precise location of the *q* barrier is known, the reactive flux of the forward reaction *k _{f}*(

*t*) can be computed from

*q*(0) −

*q*] restricts the numerator to trajectories starting at the barrier peak at time

_{b}*t*= 0, and the Heaviside function

*H*[

*q*(

*t*)] (=1 if

*q*(

*t*) >

*q*, otherwise =0) ensures that only trajectories on the product side are counted. The denominator 〈1 −

_{b}*H*(

*q*)〉 is the equilibrium probability of the reactant state. At time

*t*= 0+, one obtains

*k*=

_{f}*k*because only positive velocities can enter the product side. However, after a brief (approximately picoseconds) time period

_{TST}*τ*during which one or more barrier crossings may occur, the fate of the trajectory is determined (R or P), and the reactive flux averaged over many trajectories plateaus to a value equal to the phenomenological rate constant α (Eq. 1). Thus, the ratio of reactive flux taken at

_{mol}*t*= 0+ and

*t*=

*τ*yields the transmission coefficient κ = α/

_{mol}*k*. It is worth nothing that neither κ nor

_{TST}*k*is independent of the choice of reaction coordinate or the location of the transition state, only the product α = κ

_{TST}*k*holds that distinction, but a poor choice of reaction coordinate leads to poor statistics and possibly an erroneous estimate of α (Berne et al., 1988; Crouzy et al., 1994). Subsequent developments in reactive flux theory have improved computational speed (White et al., 2000) and broadened the application to energy landscapes with complicated topographies lacking a well-defined transition state (Dellago et al., 1999).

_{TST}Although successfully predicting the value of a rate constant by MD simulation is a notable achievement, a more informative picture of the transition process is obtained by mapping out the energy conformational landscape with phenomenological models. Similar to how experimentalists interpret voltage-clamp data using DSM models, molecular dynamicists can project the complete phase space (system plus heat bath) sampled by MD simulation onto a small number of degrees of freedom described by a low-dimensional stochastic differential equation. The generalized Langevin equation (GLE) is such an equation. Given a reaction coordinate *q* with mass *L*, the GLE describes the motion of *q* influenced by three forces: the gradient of a free energy landscape *W*(*q*), also known as the PMF, a dissipating force involving a time-dependent friction coefficient *R*(*q*,*t*), and a rapidly fluctuating force Λ(*t*). The one-dimensional GLE has the following form (Straub et al., 1987):*R*, introducing memory into the short time-scale dynamics. Λ is generally assumed to be Gaussian and is therefore completely specified by the first and second moments: 〈Λ〉 = 0 and 〈Λ(0)Λ(*t*)〉 = *kTR*(*t*). The latter equation establishes a relation between the friction-damping and fluctuating forces known as the second fluctuation-dissipation theorem (Kubo, 1966).

The task of designing MD protocols to compute *W*(*q*) and *R*(*q*,*t*) specific to electrophysiological processes has been systematically undertaken over the past two decades, much of it through the efforts of the Roux laboratory (Crouzy et al., 1994; Roux, 1997, 1999; Bernèche and Roux, 2001, 2003; Im and Roux, 2001). Recent work has sought to improve upon traditional MD force fields by incorporating polarization effects that play an important role in determining the strength of electrostatic interactions (Lopes et al., 2009; Bucher and Rothlisberger, 2010). Regardless of MD methodology, the process of coarse-graining trajectory data onto the GLE begins with the identification of a reaction coordinate *q*. The obvious (though not necessarily the most efficient for sampling purposes) choice of *q* in voltage-dependent transitions is the gating (or ionic) charge displacement, which hereafter will be understood to be synonymous with the reaction coordinate *q*. Determining the charge displacement as a function of configuration requires computing the detailed membrane potential profile in the channel with fixed (i.e., nonmobile and nonpolarizable) charges turned off (Roux, 1997). The PMF, correct to first order in *V*, is then given by *W*(*q*) = *G*(*q*) − *qV* (Roux, 1999, 2008). The zero-voltage PMF *G*(*q*) is obtained from *G*(*q*) = −*kT*ln*p*(*q*), where *p*(*q*) is the probability histogram of *q* acquired over a long MD trajectory at constant simulated temperature and pressure. Sampling statistics for the infrequently visited barrier region are improved by confining the simulation to sequential regions (umbrella windows) along the activation landscape (Torrie and Valleau, 1977; Roux, 1995). The same umbrella-sampled simulations can be used to compute *R*(*q*,*t*) from the stationary velocity autocorrelation function

Extracting the GLE parameters *W*(*q*) and *R*(*q*,*t*) from MD simulations is not a trivial undertaking, particularly if the reaction pathway (which may be multidimensional) is not well known. This is certainly true for the case of gating transitions in ion channels because available crystal structures provide only frozen “snap shots” of a particular configuration, and mechanisms of gating are poorly understood. For this reason, the early focus has been on ion permeation, where the inconvenience of dealing with an open system (grand canonical ensemble) is compensated by the existence of a relatively straightforward reaction pathway, where permeant ions have well-defined masses (Roux and Karplus, 1991). To compute rate constants from the landscape parameters *W*(*q*) and *R*(*q*,*t*), dynamic Langevin simulations can be performed (Lange and Grubmüller, 2006; Gordon et al., 2009), or the Eyring parameters ω* _{R}*, Δ

*W**, and κ can be extracted, allowing α to be calculated from Eq. 1. To execute the latter, one generally resorts to approximating the energy landscape by a piecewise sequence of parabolic basins and barriers (Fig. 1 A), treating

*R*as constant. Each segment (basin or barrier) of the landscape behaves as an oscillator satisfying the simplified GLE equation:

*t*) =

*R*(

*t*)/

*L*and the frequencies ω of local oscillations are equal to (

*LC*)

^{−1/2}. The fluctuating force ξ(

*t*) = Λ(

*t*)/

*L*satisfies 〈ξ(0)ξ(

*t*)〉 = (

*kT*/

*L*)γ(

*t*). The inverse “spring constant”

*C*is related to the curvature of the basin or barrier through

*C*

^{−1}= ε

*W″*(

*q*), where ε is 1 for (stable) basins and −1 for (unstable) barriers. Damped oscillations in the reactant and product basins (ε = 1) are exactly described by an electrical RLC circuit, where

*R*is resistance,

*L*is inductance, and

*C*is capacitance. It is straightforward to demonstrate that for a linear voltage drop over a distance

*D*, the variables

*R*,

*L*,

*C*, and

*q*are related to their mechanical equivalents (friction, mass, spring constant, and distance) through factors of

*z*/

*D*, where

*z*is the charge valence of the oscillating degree of freedom.

With Eq. 4 governing the reactant (ε* _{R}* = 1) and barrier (ε

*= −1) segments of the PMF, we easily derive ω*

_{b}*= (*

_{R}*LC*)

_{R}^{−1/2}and Δ

*W** = Δ

*G** − Δ

*q**

*V*. Obtaining the transmission coefficient κ from Eq. 4 (known as the linearized GLE because the restoring force is linear in

*q*) is more involved, but several derivations have been published (Grote and Hynes, 1980; Hänggi and Mojtabai, 1982; Pollak, 1986), including an elegant treatment based on reactive flux described by Eq. 3 (Kohen and Tannor, 2000). The satisfyingly simple result, attributed to Grote–Hynes (GH), is

*s*) is the transfer function of Eq. 4 within the barrier region and

*t*). Although evaluating λ requires specific knowledge of γ(

*t*), it is instructive to examine limiting cases. When coupling forces to the bath are weak, implying a small value of γ(

*t*) or the combination of a narrow barrier and/or small

*L*leads to a large ω

_{b}, then the rapidly accelerating

*q*at the barrier is decoupled from slower frictional forces, and the positive root of ϑ

^{−1}is λ = ω

_{b}, leading to the TST result κ = 1. (In reality the TST value represents an upper limit to the rate constant because a transition from spatial and velocity diffusion to energy diffusion occurs as

*R*is reduced to very low values, which again diminishes the value of κ [Kramers, 1940; Mel’nikov and Meshkov, 1986], although the so-called small friction regime [Fig. 1 B] is irrelevant under physiological conditions.) In contrast, the combination of a broad barrier, large

*L*, and strong friction coupling (all probable characteristics of a gating process) reduces the relative contribution of the ω

_{b}^{2}term in ϑ

^{−1}, leading to a smaller value of λ and enabling

^{−1}is easily obtained:

Kramers actually derived Eqs. 6 and 7 by implicitly assuming rapid (“memory less”) friction damping, which we can invoke by specifying γ(*t*) = 2γδ(*t*). In this case *t*)〉 = 2*kT*(γ/*L*)δ(*t*). Substituting the now time-independent γ for *k _{TST}* in Fig. 1 B). Contraction to a pure spatial diffusion (classical Brownian motion) occurs in the LF limit where the system is always at terminal velocity, enabling one to neglect the acceleration term in Eq. 8. This yields, after dividing by γ,

^{2}=

*RC*is a novel relaxation time that was implicated in the early fast component of Shaker K

^{+}channel gating currents (Sigg et al., 2003). The fluctuating force in Eq. 9 satisfies 〈

*F*(0)

*F*(

*t*)〉 = 2

*D*δ(

*t*), where the diffusion constant

*D*is related to the friction

*R*through

*D*=

*kT*/

*R*(known as the Einstein relation). Evaluating Eq. 9 for the barrier region at constant

*R*again yields

^{+}channel activation sequence (Sigg et al., 2003).

The oscillator dynamics described by Eq. 9 can be generalized to arbitrary landscape topologies *W*(*q*), which is described by the LF Smoluchowski equation, shown here in Fokker–Planck form (Van Kampen, 1992):*p*(*q*,*t*) is a probability distribution and *R*(*q*) is allowed to vary across *q*. The behavior of Eq. 11 applied to various landscape models of gating has been extensively described in Sigg et al. (1999). Operating in the spatial Smoluchowski regimen has advantages over the use of “ballistic” equations (such as the GLE), which feature velocity as well as spatial diffusion and thus require a reduced mass *L*. The absence of mass in Eqs. 9 and 11 benefits gating models because it may be difficult to define the masses of regulatory domains functioning as gating particles. Provided the energy barrier in an arbitrary {*W*(*q*), *R*(*q*)} Smoluchowski landscape is sufficiently large (>5 *kT*), then the forward transition rate constant is equal to the inverse of the mean first passage time (Van Kampen, 1992), which is easily computed from the following double integral (refer to Fig. 1 A for limits of integration):

Direct evidence for the validity of the LF assumption in gating dynamics is still lacking, although in one study of primitive gating in the gramicidin channel there was good agreement between Kramers’ LF theory and experimental rates (Crouzy et al., 1994). There is more experience with MD studies of barrier hopping of permeant and blocking ions in pores. The results from three such studies are summarized in Table 1. In each case, Langevin dynamics were found to be overdamped with κ ranging from 0.1 to 0.25. In two studies (Roux and Karplus, 1991; Tolokh et al., 2002), comparisons of κ values were made between the “gold-standard” reactive flux method (Eq. 3), the GH equation (Eq. 5), and Kramers’ intermediate-to-large friction formula (Eq. 7). There was excellent agreement between the reactive flux and GH methods, whereas Kramers’ theory underestimated κ by 20–40%, suggesting memory friction plays a nonnegligible role in ion permeation. In contrast, it is likely, based on our earlier arguments, that the larger scale of conformational changes associated with gating implies Markovian dynamics, and in the presence of extensive external (solvent, lipid) and internal (protein) damping, the Kramers–Smoluchowski condition (Eqs. 7 and 9–12) should apply. In an interesting side note, comparing Na^{+} conduction in the gramicidin pore (Roux and Karplus, 1991) with Ba^{2+} block in the KcsA channel (Rowley and Roux, 2013), the computed 10^{7}-fold difference in *k _{TST}* was primarily caused by differences in barrier height Δ

*G** because in both cases the reactant state frequency ω

*/2π was within an order of magnitude (!) of the universal rate constant*

_{R}*kT*/

*h*≈ 6 × 10

^{12}s

^{−1}.

### Comments on kinetic models

The preceding discussion sought to bridge macroscopic (DSM and TST) and microscopic (GLE, Langevin, and Smoluchowski) equations of motion by focusing on an area of common ground: the rate constant value as determined by the Eyring, reactive flux, GH, and Kramers formulas. Confidence in the validity of a phenomenological description is increased if an experimentally determined rate is consistent with the prediction of MD. However, kinetic modeling may not necessarily be the most efficient method of assigning a gating mechanism to experimental findings. To take a specific example, the mechanism by which VSDs gate the pore in voltage-dependent ion channels remains unclear. Barring direct visualization of the gating process (a distinct future possibility given advances in MD!), the researcher must infer the nature of VSD–pore coupling using indirect methods, for example by fitting ionic or gating currents to a DSM model (Vandenberg and Bezanilla, 1991; Schoppa and Sigworth, 1998b; Rothberg and Magleby, 2000) or by identifying key residues of interaction (typically highly conserved across channels) through mutation analysis or other means of functional disruption (Yang and Horn, 1995; Holmgren et al., 1998; Horn et al., 2000; Tiwari-Woodruff et al., 2000; Arcisio-Miranda et al., 2010). Model-independent methods such as mutant cycle analysis (Horovitz, 1996; Yifrach and MacKinnon, 2002; Zandany et al., 2008) can quantitatively measure coupling between residues of interacting domains such as the pore and the VSD provided that the free energies associated with the observed process are accurately measured. Mutant cycle analysis relies on conservation of energy and displacement and the ability to sum these extensive variables across a series of perturbations. Relaxation rates (eigenvalues) obtained from a multistate ion channel are not additive, although small perturbations in the activation energy of a single rate constant do add, a fact which has been used to an advantage in studying mutant cycles (Schreiber and Fersht, 1995). The problem is that, as previously stated, determining rate constants from single-channel analysis is resource intensive and typically model dependent (Colquhoun and Hawkes, 1995a; Qin, 2007). Extreme conditions are necessary to measure elementary rates directly from macroscopic currents (Sigg et al., 2003; Chakrapani and Auerbach, 2005). As a result, it is common to analyze perturbation effects through the use of standard equilibrium curves of observable quantities such as conductance (*G*) and gating charge (*Q*). This falls under the scope of thermodynamics, and although DSM models have equilibrium conditions baked into their structure, there exist simpler and more powerful methods of analyzing equilibrium curves that yield the desired energies and displacements, often in model-independent fashion. This is the subject of the next two sections.

### Thermodynamic models

In modeling electrophysiological data, thermodynamics has traditionally taken a backseat to kinetic theory, which through the use of DSM schemes enjoys widespread familiarity, but is clumsy at handling thermodynamic derivations. The thermodynamic counterpart to the Q-matrix is the partition function *Z*, a statistically weighted sum (or integral) of states from which all macroscopic quantities are derived (discussed in any statistical mechanics textbook; e.g., Hill, 1960). For a one-dimensional PMF, the partition function is*q*. A normalizing factor for *Z* is unnecessary because it merely adds a constant value to the principle quantity of interest, the channel’s chemical potential Φ, obtained from the “bridge” equation Φ = −*kT*ln*Z*. The Gibbs–Duhem equation, which describes the interdependence of intensive quantities in a thermodynamic system (Waldram, 1985), relates *d*Φ to incremental changes in system constraints (*T*, *P*, *V*, μ,…, Ω) through: *d*Φ = −〈*s*〉*dT* + 〈*v*〉*dP* − 〈*q*〉*dV* − 〈*n*〉*d*μ −…− 〈ξ〉*d*Ω, where canonical displacements 〈*s*〉, 〈*v*〉, 〈*q*〉, 〈*n*〉,…, 〈ξ〉 are written in lower case to indicate division by *N* (channels). In this way, we can derive the value of each displacement using 〈ξ〉 = −∂Φ/∂Ω = *kT*∂ln*Z*/∂Ω (an exception of sorts is entropy, which is not typically thought of as a “displacement” and possesses an extra term: 〈*s*〉 = *kT*∂ln*Z*/∂*T* + *k*ln*Z*). For example, the mean gating charge displacement is obtained from*W*(*q*) = *G*(*q*) − *qV*. We’ll examine the thermodynamics of a piecewise-harmonic two-state landscape in some detail before discussing more complex systems.

One of the many convenient properties of *Z* is the unrestricted ability to subpartition it into smaller parts, which we’ll refer to as the parsing principle. In a conformational energy landscape with large barriers, we can divide the landscape into stable (positive curvature) and unstable (negative curvature) regions. For example, the partition function of the bistable potential in Fig. 1 A may be written *Z* = *Z _{R}* +

*Z*+

_{b}*Z*, where

_{P}*X*} bounded by turning points in the curvature

*G*″(

*q*). The local potentials are second-order expansions around local extrema (

*q*,

_{X}*G*) with

_{X}*G*″ = ±

*C*

_{X}^{−1}and additionally subject to a linear voltage term −

*qV*. A large barrier minimizes the contribution of

*Z*, effectively yielding

_{b}*Z*≈

*Z*+

_{R}*Z*. After integrating and choosing (

_{P}*q*,

_{R}*G*) as the zero reference, we obtain

_{R}*W*= Δ

*G*− Δ

*qV*, where, as previously, we define Δ

*G*=

*G*−

_{P}*G*and Δ

_{R}*q*=

*q*−

_{P}*q*.

_{R}We also find it useful to propose a modified partition function that emphasizes the barrier region:*Z* by assigning a positive sign to the exponent. Performing the integral, we obtain:*W** = Δ*G** − Δ*q***V*.

The preceding exercises lay the groundwork for a rational method of discretizing a continuous energy landscape into a finite set of basin and barrier coordinates (*q _{X}*, Φ

*). These coordinates are derived from the relations:*

_{X}*q*= ±

_{X}*kT*∂ln

*Z*/∂

_{X}*V*and Φ

*= ∓*

_{X}*kT*ln

*Z*, where the upper sign applies to basins and the lower sign to barriers. Applying these expressions to Eqs. 13 and 14 yielded the coordinates of the discretized landscapes in Fig. 1 A (dashed lines). Several results ensue from this coarse-graining procedure. Revisiting kinetic theory, we see that in evaluating Eq. 12 for the bistable potential at constant

_{X}*R*, the inner integral is almost complete and approximates

*Z*at the value of

_{R}*q*where the outer integral begins to increase in value. After factoring out

*Z*, the outer integral is roughly equal to

_{R}*Z**. We therefore conclude that

_{b}*− Φ*

_{b}*and the Einstein relation was again used. Inserting Eqs. 13a and 14 into Eq. 15 and neglecting*

_{R}*V*

^{2}terms that introduce small capacity corrections to the barrier height, we see that the result is consistent with Eq. 10a obtained from Kramers’ LF theory. Also, Δ

*q** formally equals Δ

*q*

_{α}from Eq. 2a. The right side of Eq. 15 has a pleasing simplicity. The preexponential factor is the diffusion constant

*D*=

*kT*/

*R*, and the activation energy ΔΦ* is easily determined for any shape of the barrier landscape. Detailed examination of the fast component of gating current in Shaker yielded an estimate of

*D*= 41 e

_{o}

^{2}/ms, corresponding to the early stages of activation (Sigg et al., 2003).

Returning to the subject of thermodynamics, we derive the average charge displacement for the bistable potential using 〈*q*〉 = *p _{R}q_{R}* +

*p*, where

_{P}q_{P}*p*=

_{X}*Z*/

_{X}*Z*are state probabilities of the stable R and P regions. Evaluating the response to a change in voltage Δ

*V*with the aid of Eq. 13 (a and b), we obtain two components of charge movement (Fig. 1 A):

*q*=

*p*δ

_{R}*q*+

_{R}*p*δ

_{P}*q*, where the “drift” charges δ

_{P}*q*=

_{X}*C*Δ

_{X}*V*would be experimentally indistinguishable from the linear membrane capacitance if not for the fact that basin states vary in width (Sigg et al., 2003). We recall from earlier that the “fast” gating current described by δ

*q*decays with time constant

*RC*.

The second term in Eq. 16 relates to the much slower (relaxation time λ^{−1}) “hopping” or transition charge Δ*q*, which is easily distinguished experimentally from the fast component, both from a temporal standpoint (λ^{−1} >> *RC*) and also in magnitude: Δ*q* >> |δ*q _{P}* − δ

*q*| (the minus sign in the fast charge |δ

_{R}*q*− δ

_{P}*q*| is a peculiarity of the linear subtraction methods used in gating current measurements). Thus, in describing 〈

_{R}*q*〉, we typically neglect the fast charge in favor of the larger transition charge, whose equilibrium distribution is 〈

*q*〉 = [

*K*/(1 +

*K*)]Δ

*q*, where

*K*≡

*p*(∞)/

_{P}*p*(∞) =

_{R}*Z*/

_{P}*Z*is the equilibrium constant. Substituting Eq. 13 (a and b) into the above expression for

_{R}*K*, again neglecting the

*V*

^{2}terms corresponding to the fast capacity charge, we obtain

*p*(∞) = β

_{R}*p*(∞).

_{P}We note the parallel between 〈*q*〉/Δ*q* = *K*/(1 + *K*) and the venerable Boltzmann function *B* = *K _{B}*/(1 +

*K*), where

_{B}*K*= exp(ψ

_{B}*/*

_{B}*kT*) and the bias potential ψ

*= ξ*

_{B}*(Ω − Ω*

_{B}*) is a function of the force-displacement Boltzmann parameters Ω*

_{B}*and ξ*

_{B}*. The practice of Boltzmann fitting is both widely used and often abused, as shoehorning the complexity of a multidomain channel into two states requires rigorous justification (for example, Yifrach and MacKinnon, 2002). Nevertheless, Boltzmann fitting is a convenient bookkeeping technique for tracking changes to activation curves in response to systemic perturbations. Nonsymmetrical activation curves can be fitted with*

_{B}*m*Boltzmann functions, although strictly speaking the physical basis underlying such a procedure is the existence of

*m*independently gating domains described by Z = (1 +

*K*

_{1})

^{f}^{1}(1 +

*K*

_{2})

^{f}^{2}…(1 +

*K*)

_{m}*, where the*

^{fm}*f*values are population fractions. Faced with choosing between Boltzmann fitting or analysis with an

_{i}*n*-state Markov model, the Boltzmann method may in many cases be the more sensible alternative (see commentary by Shem-Ad and Yifrach, 2013). Nevertheless, a systematic approach that (a) measures thermodynamically relevant quantities, (b) uses fewer adjustable variables than a comparable kinetic model, and (c) acknowledges known or proposed structural and functional relationships between channel components is a welcomed addition to the modeling armamentarium. Such an approach (linkage analysis) has recently been applied to channel gating (Chowdhury and Chanda, 2010, 2012a, 2013; Sigg, 2013), but before discussing it, we should review ion channel thermodynamics a bit more.

One consequence of the parsing principle is the various ways *Z* may be expressed for a multistate channel. Four parsing methods (mathematically equivalent to alternative factoring schemes) applied to a linear three-state model are shown in Fig. 2. These fall into three general categories based on the method of computing the equilibrium mean for an arbitrary observable *a*:*n* and *n* − 1 states, respectively, and are useful in deriving thermodynamic relations (Conti, 1986; Sigg and Bezanilla, 1997). The relational method (Eq. 18c) differs from the other two in that the sum runs over individual transitions rather than over all states. For a linear scheme, this is no great advantage because there exist up to *n* − 1 distinct transitions in an *n*-state model, but for allosteric models possessing a small number of transitions, the dimensionality of the system can be substantially reduced. For example, the classic Monod–Wyman–Changeux (MWC) model (Monod et al., 1965) has 10 states but only two unique transitions (and one allosteric factor). A 70-state model of the BK channel (Horrigan and Aldrich, 2002) is described by three transitions and three allosteric factors. The quantity ϕ* _{K}* in relational parsing is the equilibrium curve for the transition K with range {0..

*.n*}, where

_{K}*n*is the number of identical K particles. Channel opening is typically described by a pore gate with equilibrium constant

_{K}*L*coupled to other processes (J, K). The open probability

*P*in such cases is equal to ϕ

_{o}*= ∂ln*

_{L}*Z*/∂ln

*L*, which varies from 0 to 1. It is generally straightforward to compute ϕ

*for any transition K because for DSM models,*

_{K}*Z*can always be expanded as a polynomial of

*K*.

The “particle potential” η* _{K}* in Eq. 18c, defined as the potential difference ΔΦ

*between reactant and product states of a particle, deserves special attention because it contains local thermodynamic information related to K activation. Contributions to η*

_{K}*may include particle changes in volume (Δ*

_{K}*v*) and entropy (Δ

_{K}*s*), charge movement (Δ

_{K}*q*), stoichiometries (Δ

_{K}*n*) of chemical reactions (e.g., ion, ligand, or solvent binding), and, relevant to allosterically regulated channels, energies of interaction

_{K}*W*between K and other particles J. Written explicitly, η

_{KJ}*= Δ*

_{K}*E*− Δ

_{K}*s*+ Δ

_{K}T*v*− Δ

_{K}P*q*− ΣΔ

_{K}V*n*μ

_{K}*+…+ Σ*

_{K}*W*. If one wishes to emphasize a particular extensive parameter Δξ

_{KJ}*, one can factor it out of the sum and regroup the remaining terms in the form of a corresponding half-activation value Ω*

_{K}*, yielding η*

_{K}*= −Δξ*

_{K}*(Ω − Ω*

_{K}*). For example, rewriting Eq. 17 for the equilibrium constant of a two-state system in a way that emphasizes gating charge, one obtains:*

_{K}*K*= exp(−η

*/*

_{K}*kT*) = exp[Δ

*q*(

*V*−

*V*)/

_{K}*kT*], where

*V*= Δ

_{K}*q*

^{−1}[Δ

*G*− (

*kT*/2)ln(

*C*/

_{P}*C*)]. If positive Ω activates a particle (the usual convention), then we can substitute ψ

_{R}*for −η*

_{K}*, where ψ*

_{K}*= Δξ*

_{K}*(Ω − Ω*

_{K}*) is the bias potential introduced earlier in the context of the Boltzmann function.*

_{K}The last of the four categories in Fig. 2 is the hierarchical approach, which parses *Z* in order of most to least relevant equilibrium constants (*L* is the obvious top choice given the central role of the pore). Hierarchical parsing can be seen as a relational method that leaves open the possibility of deeper categorization as new functionalities are elicited (Sigg, 2013).

### Linkage analysis

In light of the exponential dependence of an equilibrium constant *K* to its bias potential ψ* _{K}* = Δξ

*(Ω − Ω*

_{K}*), a plot of ln*

_{K}*A*versus Ω, in which

*A*=

*N*〈

*a*〉 is a marker of K activation, should yield information about Ω-linked processes relevant to K, specifically the canonical displacement Δξ

*and coupling energies*

_{K}*W*linking K to other Ω-dependent particles J. This is the basis for site-specific (local) linkage analysis developed by Jeffries Wyman to study cooperative oxygen binding in hemoglobin (Wyman, 1967) and recently adapted to ion channels (Chowdhury and Chanda, 2010, 2012a, 2013; Sigg, 2013). The parameter of interest in local linkage analysis is the Hill energy, operationally defined for an observable

_{KJ}*A*as

*W*

_{H}_{[a]}=

*kT*ln[(

*A*−

*A*)/(

_{min}*A*−

_{max}*A*)]. Defining the probability of K activation by

*P*≡

_{K}*P*= (

_{A}*A*−

*A*)/(

_{min}*A*−

_{max}*A*) and rearranging, we see that

_{min}*W*

_{H}_{[a]}=

*kT*ln[

*P*/(1 −

_{A}*P*)] is the change in negative system energy Φ from K activation: −ΔΦ

_{A}*=*

_{K}*kT*ln

*P*−

_{A}*kT*ln(1 −

*P*). Provided the particle of interest K is weakly coupled to other gating particles, a plot of

_{A}*W*

_{H}_{[a]}versus Ω is sigmoidal with linear lower and upper asymptotes that run parallel to ψ

*but are vertically separated by the total coupling energy*

_{K}*W*

_{K}_{Ω}between K and the various J particles. An example of a local variable is the pore conductance

*G*. In principle, a plot of

*W*

_{H}_{[g]}versus

*V*should yield model-independent estimates of the intrinsic pore charge Δ

*q*and the total energy of pore–VSD coupling (Sigg, 2013).

_{L}Wyman also described “global” linkage, which concerns energies of interaction between mediators of dissimilar forces (Ω = *V*, *T*, *u*_{1}, *u*_{2}, etc.; Wyman, 1967; Chowdhury and Chanda, 2013; Sigg, 2013). The parameter of interest here is the total energy of activating Ω-sensitive particles: ΔΦ_{Ω} = Δξ* _{max}*Ω

*, where Δξ*

_{M}*is total displacement related to Ω, and Ω*

_{max}*, which has been called the median (or mean, depending on interpretation) force of activation, and is defined as the value of Ω that divides the plot of 〈ξ〉 versus Ω into equal areas (the median interpretation; Wyman, 1967; Chowdhury and Chanda, 2012a), or equivalently, as the average Ω with respect to capacitance d〈ξ〉/dΩ (the mean interpretation; Di Cera and Chen, 1993; Sigg, 2013). The voltage parameter*

_{M}*V*was extensively discussed in Chowdhury and Chanda (2012a), where they demonstrated that for multistate channels marked differences may exist between

_{M}*V*and the Boltzmann-fit parameter

_{M}*V*.

_{B}The model-independent displacements (Δξ* _{K}*, Δξ

*) and energies (ΔΦ*

_{max}*, ΔΦ*

_{K}_{Ω}) obtained from local and global linkage plots are characteristic thermodynamic quantities and, when possible, should be obtained in lieu of corresponding Boltzmann parameters ξ

*and ξ*

_{B}*Ω*

_{B}*(Chowdhury and Chanda, 2012a; Bezanilla and Villalba-Galea, 2013; Sigg, 2013). The Boltzmann equation has been used as an all-purpose tool to fit any activity curve (most typically G-V or Q-V), but its validity not does extend beyond*

_{B}*n*> 2. Linkage analysis is theoretically preferable where an accurate measurement of energy or displacement is critical to the interpretation of results. A pertinent example introduced earlier is the method of double mutant cycle analysis, in which the nonadditivity of single mutant perturbation energies in double or higher dimensional mutants is considered strong evidence of interaction between mutated residues (Carter et al., 1984; Horovitz, 1987). Amino acid residues separated by 6 Å or less have been shown to interact with energies ranging from 0.5 to 7 kcal/mol (Schreiber and Fersht, 1995), with electrostatic interactions contributing to larger energies. Applied to ion channels, double mutant cycle analysis with Boltzmann fitting has been used to measure pore–VSD interactions (Yifrach and MacKinnon, 2002; Zandany et al., 2008). Chowdhury and Chanda (2010) have suggested that a local (ΔΦ

*based) linkage analysis, which they call χ analysis, be considered as a universally applicable method for quantifying mutational effects, although to date no practical application of the technique has been published, at least in part because of the technical difficulty of measuring logarithmic-scale sensitivities of local markers of activation. Chowdhury and Chanda have more recently presented a global (ΔΦ*

_{K}_{Ω}based) method of quantifying pore–VSD interactions based on the Q-V curve that was used to quantify energetic interactions in a triad of interacting residues in the Shaker K

^{+}channel (Chowdhury, S., et al. 2014. Biophysical Society 58th Annual Meeting. 3748-Pos Board B476).

Before turning to allosteric channels, a subject ideally suited to linkage analysis, we highlight a key result from Conti (1986) that has languished in relative obscurity within the ion channel literature for nearly 30 years but covers much of the thermodynamic theory described here. Starting from the equilibrium relations*p _{i}* =

*Z*/

_{i}*Z*, and ξ

*=*

_{i}*kT*∂ln

*Z*/∂Ω, Conti derived

_{i}*a*,ξ) = 〈

*a*ξ〉 − 〈

*a*〉〈ξ〉 is the covariance between any state-dependent observable

*a*and the Ω-linked displacement ξ. Eq. 19 was derived for discrete states but also covers diffusion landscapes as a consequence of the parsing principle. In the special case of

*a*= ξ, Eq. 19 simplifies to the well-known relation for capacitance: kT∂〈ξ〉/∂Ω = 〈Δξ

^{2}〉. The Conti relation is insensitive to the number of channels

*N*because replacing 〈

*a*〉 with

*A*=

*N*〈

*a*〉 leaves Eq. 19 unchanged. We can equate the open probability

*P*with the unit-normalized G-V curve = (

_{o}*G*−

*G*)/(

_{min}*G*−

_{max}*G*), where

_{min}*G*=

*N*〈

*g*〉 and

*G*is a baseline conductance. Applying Eq. 19 to

_{min}*P*with Ω =

_{o}*V*, we obtain an expression for the “activation charge” 〈

*q*〉 ≡

_{a}*kT*dln(

*G*−

*G*)/d

_{min}*V*= 〈

*q*〉 − 〈

_{o}*q*〉, where 〈

*q*〉 = 〈

_{o}*P*〉/

_{o}q*P*is the mean charge displacement among open states (Sigg and Bezanilla, 1997). Similarly, defining the closed probability

_{o}*P*= 1 −

_{c}*P*, we obtain the complementary quantity 〈

_{o}*q*〉 ≡ −

_{d}*kT*dln(

*G*−

_{max}*G*)/d

*V*= 〈

*q*〉 − 〈

*q*〉, where 〈

_{c}*q*〉 = 〈

_{c}*P*〉/

_{c}q*P*is the closed-state charge displacement, and we also have 〈

_{c}*q*〉 =

*P*〈

_{c}*q*〉 +

_{c}*P*〈

_{o}*q*〉. Finally, we recognize that for positive voltage-activating channels the slope of the conductance Hill plot

_{o}*m*= d

*W*

_{H}_{[g]}/d

*V*≡ 〈

*q*〉 + 〈

_{a}*q*〉 = 〈

_{d}*q*〉 − 〈

_{o}*q*〉 asymptotically approaches Δ

_{c}*q*for extreme negative (

_{L}*m*

_{(−)}≈ 〈

*q*〉) and positive (

_{a}*m*

_{(+)}≈ 〈

*q*〉) voltages because under these limiting conditions, the voltage sensor is locked into either the resting or activated position and the pore can gate quasi-independently. Channels for which

_{d}*m*

_{(−)}=

*m*

_{(+)}measured under experimental conditions may be said to exhibit “weak” pore–VSD coupling, and gating schemes can either conform to this condition or violate it through “strong” or “obligatory” coupling, as discussed presently.

### Allosteric models of gating

Some 10 years after Hodgkin and Huxley’s groundbreaking work on the action potential, a different, nonkinetic approach yielding the well-known MWC allosteric model was applied to the problem of oxygen binding in hemoglobin. The MWC scheme proposed that the O_{2} affinities of the *n* = 4 hemoglobin subunits increased simultaneously with the flip of a two-state L parameter (Monod et al., 1965). The ligand-gated ion channel community was quick to catch on, with publication two years later of an *n* = 2 MWC model for the acetylcholine channel (Karlin, 1967). The physical basis of L in both hemoglobin and ligand-gated channels is thought to lie with symmetry-preserving residue interactions at the interfaces between subunits (Changeux and Edelstein, 2005). In voltage-dependent ion channels, a case can be made equating the pore with the “L particle” because pore opening is a binary process *n*-fold connected to neighboring VSDs. It wasn’t until the 1990s that allosteric models of voltage-gated channels were considered seriously, culminating in two influential papers from the Aldrich laboratory. The first (Ledwell and Aldrich, 1999) studied a late opening transition in the voltage-sensitive Shaker K^{+} channel by mutating three noncharged residues of the mid-to-distal portion of the VSD S4 segment that interact with the pore domain (Pathak et al., 2007). The authors found that this so-called ILT mutant (Smith-Maxwell et al., 1998) led to single-exponential time courses in the ionic current and a far right-shifted G-V curve, consistent with a rate-limiting opening transition. They interpreted their findings in the context of traditional obligatory gating models of Shaker activation (Bezanilla et al., 1994; Zagotta et al., 1994; Schoppa and Sigworth, 1998b) in which pore opening occurs only at the end of the activation sequence, as originally prescribed by the HH scheme. The apparent absence of more than one open state in obligatory gating is at odds with the notion of weak allosterism, but a kind of strong coupling is inferred by the requirement that all four VSDs must activate before pore opening. Support for obligatory gating in Shaker and other homologous voltage-gated K^{+} channels came from “limiting slope” experiments (Almers, 1978; Sigg and Bezanilla, 1997), in which the activation charge 〈*q _{a}*〉 experimentally approaches the value of total activation charge Δ

*q*at very low

_{max}*P*∼10

_{o}^{−7}(Islas and Sigworth, 1999). At such low

*P*, the gating charge versus voltage curve (Q-V) is negatively saturated (〈

_{o}*q*〉 ≈ 0) and 〈

*q*〉 =

_{a}*m*

_{(−)}= Δ

*q*). Evidently, the open state charge 〈

_{max}*q*〉 remains maximal in Shaker even at very low open probabilities, as demonstrated by 〈

_{o}*q*〉 ≡ 〈

_{o}*q*〉 + 〈

_{a}*q*〉 ≈ Δ

*q*. This rules out the existence of intermediate-charged open states, for which 〈

_{max}*q*〉 < Δ

_{o}*q*. Obligatory gating in Shaker is consistent with single-exponential open time distributions in singles analysis (Hoshi et al., 1994; Schoppa and Sigworth, 1998a) and the observation that just one charged voltage sensor is sufficient to gate the pore shut (Gagnon and Bezanilla, 2009). However, if one takes the position that obligatory gating represents the strong coupling limit of an underlying allosteric scheme (a position which is admittedly at odds with current conventional wisdom), then one obtains a lower estimate of the pore–VSD coupling energy

_{max}*W*required to approach the limiting slope

*m*

_{(−)}= Δ

*q*for

_{max}*P*= 10

_{o}^{−7}by assuming that the corresponding point in

*W*

_{H}_{[g]}occurs halfway between lower and upper asymptotes (corresponding to the steepest portion of the Hill plot), yielding

*W*= 2

*kT*ln(10

^{−7}) ≈ 32

*kT*, or ∼8

*kT*per subunit. An allosteric model of Shaker consistent with the ILT experiments that nearly achieves this coupling requirement is presented later.

The second of the aforementioned Aldrich papers (Horrigan and Aldrich, 2002), which focused on the large conductance (BK) Ca^{2+}- and voltage-activated K^{+} channel, is arguably the most complete characterization of allosterism in a voltage-gated cation channel to date. The energies of interaction between the VSD (J particle), the pore (L particle), and the Ca^{2+}-binding domain (K particle) obtained in this landmark study ranged from 0.5 to 1.9 kcal/mol (0.9 to 3.2 *kT*) per subunit, with the upper limit in energy corresponding to J-L coupling. The negative limiting slope *m*_{(−)} in BK for *P _{o}* < 10

^{−3}(high [Ca

^{2+}]) and 10

^{−7}(low [Ca

^{2+}]) yielded Δ

*q*= 0.3 e

_{L}_{o}. Weak allosterism has been implicated in gating models of other channels, with derived interaction energies ranging from 1.0 to 1.2 kcal/mol (1.7 to 2.0

*kT*) in the voltage- and temperature-gated TRPM8 channel (Brauchi et al., 2004), 0.9 to 1.8 kcal/mol (1.6 to 3.2

*kT*) in voltage- and cyclic nucleotide-gated HCN channels (Altomare et al., 2001; Chen et al., 2007; Ryu and Yellen, 2012), 0.4 to 1.0 kcal/mol (0.7 to 1.8

*kT*) in the voltage-gated KCNQ1 channel (Osteen et al., 2012), 2.9 kcal/mol (5.0

*kT*) in a sodium channel (Arcisio-Miranda et al., 2010), and 1.6 kcal/mol (2.7

*kT*) in a calcium channel (Marks and Jones, 1992).

The basic unit of interaction in an allosteric scheme is the pairwise coupling between two gating particles that we’ll refer to as a “cooperon.” Consider an isolated cooperon between J and L particles whose Gibbs energy landscapes are shown in Fig. 3 A. The corresponding single particle partition functions are *Z _{J}* = 1 +

*J*and

_{o}*Z*= 1 +

_{L}*L*, where

_{o}*J*and

_{o}*L*represent native equilibrium constants. Interactions between the activation states of J and L, if not too strong or distorting, generate a two-dimensional landscape of four basins arranged in a cycle (Sigg and Bezanilla, 2003). The partition function in each basin state can be constructed from the particle states (

_{o}*j*,

*l*) and their strengths of interaction

*D*. Thus, the cooperon partition function is

_{jl}*Z*=

*D*

_{00}+

*J*

_{o}D_{10}+

*L*

_{o}D_{01}+

*J*

_{o}L_{o}D_{11}. The interaction free energies and charges are given by

*W*= −

_{jl}*kT*ln

*D*and δ

_{jl}*q*=

_{jl}*kT*∂ln

*D*/∂

_{jl}*V*, with the latter usually set to zero. Dividing

*Z*by

*D*

_{00}and reassigning terms yields the more compact expression

*Z*= 1 +

*J*+

*L*+

*JLD*, where

*J*=

*J*

_{o}D_{10}/

*D*

_{00},

*L*=

*L*

_{o}D_{01}/

*D*

_{00}, and

*D*=

*D*

_{00}

*D*

_{11}/

*D*

_{10}

*D*

_{01}are renormalized variables derived from native equilibrium constants and raw coupling factors (Chowdhury and Chanda, 2010). There are no restrictions on manipulations of this type because, as noted earlier, normalizing factors do not change energy relations between states. Because the global coupling factor

*D*targets the double product state A-O (Fig. 3), we refer to the renormalized scheme as “product coupled.” The fact that

*D*is a composite of all four native interactions implies that an apparent absence of interaction between two particles (

*D*= 1) does not necessarily imply complete particle independence, but may be a consequence of balanced energies, specifically

*W*≡

*W*

_{00}+

*W*

_{11}−

*W*

_{10}−

*W*

_{01}= 0. We also consider the possibility of excess charge δ

*q*=

*kT*∂ln

*D*/∂

*V*, which renders

*D*voltage dependent provided that δ

*q*= δ

*q*

_{00}+ δ

*q*

_{11}− δ

*q*

_{10}− δ

*q*

_{01}≠ 0. Evidence has been offered in support of a small (0.2 e

_{o}) δ

*q*in internal coupling between BK voltage sensor particles (Pantazis et al., 2010).

In most voltage-gated ion channels, *P _{o}* increases with voltage. A physical mechanism consistent with this property is when voltage sensors assist the opening of a reluctant pore gate through positive interaction (negative

*W*) in the A-O state. This “pull-open” mechanism is illustrated on the right side of Fig. 3 B; it satisfies the product-coupled partition function

*Z*= 1 +

*J*+

*L*+

*JLD*, with

*D*> 1. However, one can also imagine a “push-closed” mechanism (Fig. 3 B, left) in which the pore opens intrinsically, but only after activation of the J particle relieves a steric hindrance

*D*′ in the R-O state, which prevents pore opening while J is at rest (Chowdhury and Chanda, 2012b; Horrigan, 2012). The push-closed mechanism is described by

*Z*= 1 +

*J*+

*L*′

*D*′ +

*JL*′ with

*D*′ < 1. It is entirely equivalent to the pull-open mechanism provided we set

*D*′ = 1/

*D*and

*L*′ =

*LD*. The conductance Hill plots of the two mechanisms coincide if their respective L particle biases ψ

*and ψ*

_{L}*′ are vertically separated by the coupling energy W = −*

_{L}*kT*ln

*D*(Fig. 3 B). In either case, if we allow

*L*to become very small and

*D*very large, then

*Z*≈ 1 +

*J*+

*JL*′ (Arcisio-Miranda et al., 2010), which describes the three-state obligatory model of Fig. 2 and represents the limiting case of L-J coupling in Fig. 3 B (dotted red line). The extension to a multi-cooperon scheme consisting of a central pore particle L radially coupled to

*n*J particles is straightforward, where

*Z*= (1 +

*J*)

*+*

^{n}*L*(1 +

*JD*)

*for the pull-open mechanism,*

^{n}*Z*= (1 +

*J*)

*+*

^{n}*L*′(

*D*′ +

*J*)

*with*

^{n}*L*′ =

*LD*for the push-closed mechanism, and

^{n}*Z*= (1 +

*J*)

*+*

^{n}*L*′

*J*for the (

^{n}*n*+ 2)-state obligatory model.

If, practically speaking, one could eliminate the source of cooperativity *D* (or *D*′) from a product-coupled scheme, then the ensuing position of *V _{L}* (or

*V*′) on the

_{L}*V*axis should favor a particular gating mechanism (pull-open or push-closed). It is in this context that an alternative interpretation of the ILT experiments is proposed, one that regards Shaker gating as the manifestation of strong coupling rather than following a strictly obligatory scheme. We entertain the possibility that the ILT mutation completely uncouples the pore from the voltage sensor, reducing the value of

*D*to the neutral value 1.0. In fact, some residual pore–VSD coupling may exist in the ILT mutant (Smith-Maxwell et al., 1998), but we’ll treat this as negligible. The observed G-V shift to the right in ILT experiments should therefore favor the pull-open hypothesis. The difficulty with this interpretation relates to the change in the Q-V curve in response to ILT, which undergoes an opposite shift to the left (Ledwell and Aldrich, 1999). It is generally the case that changes to a single coupling factor or equilibrium constant is inconsistent with opposing activation shifts (Chowdhury and Chanda, 2012b). Similar examples of diverging G-V and Q-V curves have since been observed in a handful of Shaker mutants, mostly involving residues found in the interaction region between the VSD and the pore (segments S5 and S6 in the pore and the S4–S5 linker; Soler-Llavina et al., 2006; Haddad and Blunck, 2011). Interestingly, neutralizing salt bridge–forming charged residues in the S2 and S3 segments of the VSD (Papazian et al., 1995; Tiwari-Woodruff et al., 1997) have also generated opposing shifts (Seoh et al., 1996).

Muroi et al. (2010) provided a plausible mechanism to explain the divergence of G-V and Q-V curves upon observing the phenomena in sodium channel mutants targeting the domain III VSD–pore linkage region. They showed that if a mutation were to simultaneously affect separate coupling factors in the L-J cooperon, specifically if one factor targets the doubly “reactant” state (R-C) and the other the “product” state (A-O), a scheme referred to here as “doubly coupled,” then relieving both sources of coupling could generate opposing shifts in ϕ* _{L}* and ϕ

*and also therefore in the G-V and Q-V curves. The concept of interactions favoring “agreement” states such as R-C and A-O is a common one in statistical mechanics. It explains the quaternary symmetry requirements of the MWC model (Monod et al., 1965) and underlies the Lenz–Ising model of lattice interactions (Niss, 2005).*

_{J}Incorporating double coupling into an *n* = 4 radial cooperon scheme (Fig. 4) yields the partition function*L* and *J* are the pore and voltage sensor equilibrium constants, and the coupling factors associated with R-C and A-O configurations are *C* = exp(−*W _{C}*/

*kT*) and

*D*= exp(−

*W*/

_{D}*kT*), respectively. Assuming in Shaker a total charge per channel Δ

*q*= 13 e

_{max}_{o}(Aggarwal and MacKinnon, 1996; Noceti et al., 1996; Seoh et al., 1996; Islas and Sigworth, 1999), Eq. 20 adequately describes the changes made by the ILT mutations to the G-V and Q-V curves using just five adjustable parameters (Δ

*q*,

_{L}*V*,

_{L}*V*,

_{J}*W*, and

_{C}*W*), all related to displacements and energies obtained from linkage analysis (Fig. 4, B and C).

_{D}Although the *n* = 4 doubly coupled scheme is allosteric, with half of its 32 states conducting, the ILT fits account for obligatory-like gating in Shaker (Fig. 4 C) by distributing the total energy of interaction *W _{tot}* = 4|

*W*+

_{C}*W*| ≈ 26

_{D}*kT*equally among the four subunits, with 70% of

*W*allotted to R-C interactions (4|

_{tot}*W*| ≈ 18

_{C}*kT*) and 30% to A-O interactions (4|

*W*| ≈ 8

_{D}*kT*). This value of

*W*is roughly 80% of the earlier back-of-the-envelope estimate of 32

_{tot}*kT*for the minimum VSD–pore coupling energy required to emulate obligatory gating in Shaker, leaving open the possibility that a more completely uncoupled mutant (see Haddad and Blunck, 2011) might yield a larger value of

*W*through greater separation of G-V and Q-V curves. We note that mathematically there is no difference between stabilizing the double reactant and product configurations and destabilizing intermediate states. Therefore, ILT or similarly acting mutations could conceivably act by weakening favorable interactions in agreement states R-C and A-O as suggested above or by relieving “strain” in intermediate states R-A and C-O (Pathak et al., 2005; Haddad and Blunck, 2011).

_{tot}One could argue that, by factoring out the quantity *C*^{4} from Eq. 20 to obtain *Z* = (1 + *J*′)^{4} + *L*′(1 + *J*′*CD*)^{4}, with effective coupling constant *CD* and redefined equilibrium constants *J*′ = *J*/*C* and *L*′ = *L*/*C*^{4}, one simply reverts back to a product-coupled scheme of the pull-open or reformulated push-closed type. This is true, but the above function differs from those of the earlier schemes in that both *J*′ and *L*′ are decreased in value with increasing *C* (recall that *J* in the pull-open or push-closed schemes was unchanged in the limit of large coupling). Hence, if *C* → ∞, the doubly coupled scheme reduces to a two-state model described by *Z* = 1 + *LJ*^{4}(*D*/*C*)^{4}, provided we also have a large *D* to maintain energy balance (note the absence of opening bias if *D* = *C*). A limiting doubly coupled scheme of this type might explain the binary activation of the four-subunit pore gate, which prefers all S6 helical bundle segments to be in the closed or open state, thereby satisfying quaternary symmetry required by MWC theory (Chapman and VanDongen, 2005).

A mechanistic interpretation compatible with the doubly coupled scheme is if disruption of packing interactions between the pore and VSD is necessary for intermediate gating transitions to occur, a conjecture supported by the findings of long MD trajectories of K^{+} channel activation (Jensen et al., 2012). Although the doubly coupled scheme is overly simplistic, neglecting for example the multistep activation of the VSD (Zagotta et al., 1994; Schoppa and Sigworth, 1998b; Gamal El-Din et al., 2010; Delemotte et al., 2011; Haddad and Blunck, 2011; Yarov-Yarovoy et al., 2012), it adopts a thermodynamic viewpoint in which the parameters of interest consist solely of energies and displacements. This lines up perfectly with the way we often think of mechanisms of action, which focus on mechanical details such as which residues in domain X interact with those in domain Y and how these interactions are modified when there is a change in conformation. Kinetic modeling provides additional information about the time course of events, but this may not be of central importance in determining structure–function relationships, and the added complexity of a kinetic model can obscure interpretation. The double coupling framework is easily expanded to include multistep VSD transitions and could be used to explore other features of gating, such as intersubunit pore domain interactions and channel inactivation. Adding new gating particles and transitions to a multimeric allosteric model can result in an unwieldy proliferation of kinetic states, although equilibrium curves can still easily be calculated using the partition function formulism. Of course, a successful model must be able to predict gating and ionic currents, and for this purpose it is necessary to adopt a systematic approach to kinetic analysis of allosteric models, as discussed next.

### Kinetics in allosteric models

Constructing a kinetic model of allosterism can be a daunting task. Thermodynamics offers only incomplete guidance on constraints placed by allosteric interactions on the perturbation of rate constants, and the number of kinetic states rises geometrically with cooperon number. A reasonable solution to the first problem is the concept of linear free energy relationships (LFERs; Leffler, 1953; Hammond, 1955; Grosman et al., 2000; Azaria et al., 2010; Edelstein and Changeux, 2010). The basic idea is that if, structurally speaking, the transition state (b) is intermediate between the reactant (R) and product (P) states, then the perturbation of state potentials Φ* _{X}* by an arbitrary force Ω should satisfy the weighted average: ∂Φ

*/∂Ω =*

_{b}*x*∂Φ

*/∂Ω + (1 −*

_{P}*x*)∂Φ

*/∂Ω, where the weighting factor*

_{R}*x*(LFER parameter) is a number between 0 and 1. Rearranging terms, we have: ∂ΔΦ*/∂Ω =

*x*∂ΔΦ/∂Ω, where ΔΦ* = Φ

*− Φ*

_{b}*and ΔΦ = Φ*

_{R}*− Φ*

_{P}*. We obtained earlier (Eq. 15) α =*

_{R}*D*exp(−ΔΦ*/

*kT*), from which we derive, assuming Φ is linear in Ω, α(Ω)/α(0) ≡ exp[−(∂ΔΦ*/∂Ω)Ω/

*kT*] = exp[−

*x*(∂ΔΦ/∂Ω)Ω/

*kT*]. For the familiar case of Ω =

*V*and substituting ∂ΔΦ/∂

*V*= −Δ

*q*, we obtain α = α

_{o}exp(

*x*Δ

*qV*/

*kT*), where, comparing with Eq. 2a, we establish that Δ

*q*(=Δ

_{a}*q**) =

*x*Δ

*q*. Extending the above argument to both forward and backward rates, we can write

*= exp(Δ*

_{V}*qV*/

*kT*) is the voltage perturbation factor. The ratio α/β yields the usual expression for the equilibrium constant K =

*K*exp(Δ

_{o}*qV*/

*kT*), where

*K*= α

_{o}_{o}/β

_{o}. We can introduce additional factors for other applied forces (Γ

*, Γ*

_{T}*, Γ*

_{P}_{μ},…) as well as for allosteric interactions (

*j*,

*l*,… are the number of activated J, L,… particles. Generally speaking, the LFER parameter

*x*of a particular perturbation ∂Ω depends on the position of the transition peak along the coordinate ξ = −

*kT*∂Φ/∂Ω. However, we can make the zero-order assumption that

*x*has the same value for any perturbation and further specify a transition-specific relaxation rate ν = α

_{o}

^{(1−x)}β

_{o}

^{x}, allowing us to economically write: α = ν

*K*and β = ν

^{x}*K*, where

^{x−1}*= −*

_{K}*kT*ln

*K*, we have

*K*, α, and β, the resulting relations represent the basis for a phenomenological rate description. Combinatorial factors (σ

*= 1…*

_{K}*n*) must be included as part of a global kinetic scheme if there are

*n*identical K particles. A heuristic model of the BK channel comprising

*m*= 4 particle species arranged in fourfold symmetry around the pore for a total of 17 particles and 7,752 configurational states (twice the number of ways

*n*= 4 subunits can occupy 2

*= 16 subunit states) was recently evaluated in this manner (Sigg, 2013), although the large state space made kinetic analysis computationally intensive, requiring Monte Carlo methods.*

^{m}A shortcut bypassing the rigorous Q-matrix procedure of kinetic analysis would be of considerable benefit when dealing with large-scale allosteric models. The earlier trick of computing equilibrium curves by summing over elementary transitions using the relational form of the partition function (Eq. 18c) suggests one might construct an analogous set of coupled HH-like kinetic equations, one for each gating particle, where the discrete coupling numbers *j*, *l*,… in Eq. 21 are replaced with the time-varying ϕ* _{J}*, ϕ

*,… in the form of a mean-field approximation. Unfortunately, the Q-matrix of allosteric models is not so easily coarse grained. Although the mean-field procedure satisfies equilibrium requirements, it distorts the time course of relaxation even for the simple cooperon scheme. However, if one particle is rate limiting, a successful approach analogous to the Born–Oppenheimer strategy from quantum chemistry can be used. Consider an*

_{L}*n*-meric channel with a rate-limiting pore L relative to other regulatory particles (ν

*<< ν*

_{L}*, ν*

_{J}*,…). The partition function can be written as*

_{K}*Z*=

*Z*+

_{c}^{n}*LZ*, where the closed-state component

_{o}^{n}*Z*=

_{c}*f*(

*J*,

*K*,…) is a polynomial function of subunit equilibrium constants

*J*,

*K*,…. The same polynomial describes the open state partition function

*Z*=

_{o}*f*(

*JD*,

*KC*,…) except that the equilibrium constants are multiplied by allosteric factors

*D*,

*C*,…. The pore then activates in the manner of an HH particle, characterized by

_{L}_{∞}= ∂ln

*Z*/∂ln

*L*and λ

*= α*

_{L}*+ β*

_{L}*. The rate constants α*

_{L}*and β*

_{L}*are Boltzmann-weighted sums of elementary pore rates across {J, K,…} substates, each individually constructed in the manner detailed in the previous paragraph. Evaluating the two sums yields the compact expression*

_{L}*Z*=

_{x}*f*(

*JD*,

^{x}*KC*,…) is a sliding partition function which varies from

^{x}*Z*to

_{c}*Z*in the range

_{o}*x*= {0…1}. The fast gating particles are “slaved” to L and therefore relax with rate λ

*after an initial rapid decay. An approach mathematically analogous to Eq. 22 was used to compute ionic current relaxation rates in the BK channel, whose pore kinetics are roughly 10 times slower than those of its regulatory domains (Cox et al., 1997; Horrigan and Aldrich, 2002).*

_{L}### Conclusion

The history of modeling ion channels is filled with many successes and some dead ends. Even the dead ends are useful in that they test the limits of our understanding. This review has attempted to bridge microscopic and macroscopic views of ion channel dynamics through an exploration of common ground approaches to the phenomenology of gating. Although the lion’s share of attention was directed toward voltage-dependent channels, the methods discussed here are easily generalized to other stimuli and are applicable to any channel and indeed any macromolecule. Looking forward, close cooperation between the disciplines of MD and electrophysiology will likely play an important role in unraveling the molecular mechanism of gating. Along these lines, a consensus view of the activation pathway for voltage-gated ion channels based on decades of experimentation and simulation modeling has recently been published (Vargas et al., 2012). A logical next step would be to map out potentials of mean force for the gating landscape under physiological conditions. This would require a large effort (some might call it a pipe dream) but offers a substantial payoff: the validation, based on a channel’s molecular structure, of specific gating models used with confidence to interpret experimental data, with the ultimate aim of improved drug development and alleviation of disease. A combined approach has recently been applied (Ostmeyer et al., 2013) to the question of why recovery from C-type inactivation in the K^{+} pore is so slow (∼5–20 s). The answer? Three water molecules trapped behind the selectivity filter of each subunit bolster the pinched conformation of the inactivated pore. PMF calculations in this study demonstrated a large energy barrier (∼25 kcal/mol) to pore conduction unless all (or most) of the 12 waters were released into solution. With the benefit of hindsight and using thermodynamic arguments, one might have guessed the role of inactivating waters based on earlier work that studied the pressure dependence of slow inactivation in Shaker (Meyer and Heinemann, 1997), but molecular simulation combined with kinetic modeling was necessary to illustrate the precise mechanism in nearly cinematic detail. We should expect more such exciting developments arising from the interface between molecular and macroscopic dynamics. Whether an accurate phenomenology of gating dynamics follows the script suggested by the contents of this review, or whether new modeling paradigms come into play, based in part on advances in computer simulations but informed by experiment, should become clearer over the next decade or so as more light is shed on this fascinating and important issue.

## Acknowledgments

Helpful comments on the manuscript by Dr. Baron Chanda and Dr. Richard Aldrich are gratefully acknowledged.

The author declares no competing financial interests.

Kenton J. Swartz served as editor.

## Footnotes

- Abbreviations used in this paper:
- DSM
- discrete-state Markov
- GH
- Grote–Hynes
- GLE
- generalized Langevin equation
- HH
- Hodgkin–Huxley
- LF
- large friction
- LFER
- linear free energy relationship
- MD
- molecular dynamics
- MWC
- Monod–Wyman–Changeux
- PMF
- potential of mean force
- TST
- transition state theory
- VSD
- voltage-sensing domain

- Submitted: 30 October 2013
- Accepted: 21 May 2014

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).