Electrophysiological studies have established that the permeation of Ba2+ ions through the KcsA K+-channel is impeded by the presence of K+ ions in the external solution, while no effect is observed for external Na+ ions. This Ba2+ “lock-in” effect suggests that at least one of the external binding sites of the KcsA channel is thermodynamically selective for K+. We used molecular dynamics simulations to interpret these lock-in experiments in the context of the crystallographic structure of KcsA. Assuming that the Ba2+ is bound in site S2 in the dominant blocked state, we examine the conditions that could impede its translocation and cause the observed “lock-in” effect. Although the binding of a K+ ion to site S1 when site S2 is occupied by Ba2+ is prohibitively high in energy (>10 kcal/mol), binding to site S0 appears to be more plausible (ΔG > 4 kcal/mol). The 2D potential of mean force (PMF) for the simultaneous translocation of Ba2+ from site S2 to site S1 and of a K+ ion on the extracellular side shows a barrier that is consistent with the concept of external lock-in. The barrier opposing the movement of Ba2+ is very high when a cation is in site S0, and considerably smaller when the site is unoccupied. Furthermore, free energy perturbation calculations show that site S0 is selective for K+ by 1.8 kcal/mol when S2 is occupied by Ba2+. However, the same site S0 is nonselective when site S2 is occupied by K+, which shows that the presence of Ba2+ affects the selectivity of the pore. A theoretical framework within classical rate theory is presented to incorporate the concentration dependence of the external ions on the lock-in effect.
K+ channels are a broad family of membrane proteins that are present in almost every cell. Their high selectivity is one of the most remarkable aspects of cellular physiology. While remaining highly selective for K+ ions over Na+ ions, these channels allow K+ ions to cross a cellular membrane through a passive mechanism at a rate that nearly matches bulk diffusion. Conduction occurs through a transmembrane domain with tetrameric architecture. The subunits of this tetrameter contain a highly conserved sequence of amino acids, TTVGYGD, which is essential for K+ selectivity (Heginbotham et al., 1994). The high-resolution crystal structure of KcsA, a pH-gated K+ channel from the bacteria Streptomyces lividans, revealed that the selectivity filter provides a narrow pore for single-file ion conduction comprised of a series of ion binding sites where the ion is coordinated by backbone carbonyl oxygens (Doyle et al., 1998). This structure is illustrated in Fig. 1. The robust selectivity achieved by such a structurally simple motif is even more intriguing considering the sizable structural fluctuations that proteins undergo at physiological temperature (Allen et al., 2004). This is especially remarkable given that the sizes of K+ and Na+ ions are very similar; the ionic radius of K+ is only ∼0.4 Å larger than Na+ (Shannon, 1976).
Under physiological conditions, the open-conductive state of the channel must allow K+ ions to move outwards from the intracellular solution to the extracellular solution without allowing Na+ ions to enter from the extracellular solution. In effect, this means that the narrow filter must select for K+ ions in its open-conductive conformation, without undergoing a long-timescale conformational change.
Since the report of the high-resolution x-ray structure, KcsA has become a prototypical model for experimental as well as computational studies of the structure and activity of K+ channels. A series of perspectives on ion channel selectivity was recently published in this journal (Andersen, 2011). Many computational studies have attempted to explain the structure, dynamics, and selectivity of the filter (Sansom et al., 2002; Shrivastava et al., 2002; Bernèche and Roux, 2003; Domene et al., 2008; Domene and Furini, 2009; Furini and Domene, 2009). In particular, there have been many efforts to elucidate the origin of selective ion conduction using molecular dynamics simulations (Allen et al., 2000; Aqvist and Luzhkov, 2000; Shrivastava et al., 2002; Noskov et al., 2004; Roux, 2005; Egwolf and Roux, 2010; Jensen et al., 2010; Kim and Allen, 2011). These studies have led to a range of interpretations of ion channel selectivity, with Roux and coworkers concluding that K+ selectivity is enforced by a higher free energy for Na+ ions binding in site S2, whereas Kim and Allen (2011) concluded that no binding site is significantly selective for K+. In part, the myriad of theories on the origin of the K+ channel selectivity result from difficulty in relating these simulations to the experimentally observed activity of these channels.
One of the most direct measures of K+ selectivity is provided by blockade experiments in the presence of Ba2+ ions (Hagiwara et al., 1978; Armstrong and Taylor, 1980; Neyton and Miller, 1988; Vergara et al., 1999; Piasta et al., 2011). As the radius of Ba2+ is very similar to that of K+ (1.38 Å and 1.35 Å, respectively), both of these cations can fit inside the binding sites and permeate through the selectivity filter of a K+ channel. However, the divalent Ba2+ ions bind more strongly to the selectivity filter than K+ ions and, thus, permeate through the channel at a much slower rate. The binding of Ba2+ to the filter prevents the passage of other ions, resulting in blockades of the channel current on the millisecond timescale that are observed in single-channel recordings. Recently, such Ba2+ blockade experiments were performed for the first time with the KcsA channel (Piasta et al., 2011). Such experiments are expected to be particularly informative in this case because of the availability of high-resolution x-ray structure of this channel. When a 30−60-µM concentration of Ba2+ ions was introduced in the “internal” solution, the single channel current underwent frequent interruptions due to Ba2+ binding to the selectivity filter. The block is relieved when the Ba2+ ion exits the selectivity filter, either by returning to the intracellular side or proceeding to the external solution. Block times were typically of the order of 1 ms in the absence of external K+, but increased to ∼100−1,000 ms with an external K+ concentration of 5 mM. Even ∼5-µM external concentrations of K+ increase the duration of the Ba2+ blockade lifetimes, which suggests the existence of an outer binding site with a high affinity for K+. In contrast, experiments with Na+ ions in the external solution showed no lock-in effect at all at concentrations up to 0.2 M. This dependence of the block time on the concentration of extracellular ion is interpreted as evidence of the so-called K+ lock-in effect, where the binding of K+ in a site external to the Ba2+ prolongs the blocking time by impeding its translocation forward toward the extracellular side. Because the blocking times are very long, these experiments can be interpreted by assuming that binding of a cation to a site located external relative to the bound Ba2+ is in near thermodynamic equilibrium with the extracellular bulk solution. Using a kinetic model based on this view, it was deduced that the external lock-in site must be thermodynamically selective for K+ over Na+ by at least 7 kcal/mol.
In many respects, these Ba2+ blockade experiments are highly complementary to detailed free energy computations using molecular dynamics simulations. As the permeation of Ba2+ occurs on a time scale that is many orders of magnitude larger than the rate of Na+ or K+ association and dissociation from the filter, the alkali ions binding in the external sites of the filter can be safely understood as being in equilibrium with the external solution. In this way, the selectivity of the channel does not need to be discussed in terms of the nonequilibrium diffusive process of conduction, but rather it is framed in terms of the equilibrium thermodynamic binding affinity of the external sites to an alkali ion and the quasi-equilibrium process of rare, activated transitions of Ba2+ between sites. The Ba2+ blockade experiments form the basis of a thermodynamic view of ion selectivity, where differences in the binding affinity of K+ or Na+ for specific sites along the narrow pore is the mechanism by which these channels achieve selective ion permeation. In this paper, we report umbrella sampling and free energy perturbation (FEP) molecular dynamics simulations designed to model the Ba2+ blockade experiments of Piasta et al. (2011) in an attempt to reconcile the experimental observations with the results of simulation.
MATERIALS AND METHODS
Models and computational parameters
All simulations reported here were performed using the program CHARMM, version c36b2 (Brooks et al., 2009). The atomic simulation systems were based on the systems that were used in previous simulations by our group (Bernèche and Roux, 2000; Egwolf and Roux, 2010). The tetrameric crystallographic structure of the KcsA channel (Protein Data Bank [PDB] accession no. 1K4C) reported by Zhou et al. (2001) is embedded in a 70 Å × 70 Å bilayer comprised of 112 dipalmitoyl-phosphatidylcholine (POPC) lipids. The pore axis of the channel was aligned along the z axis of the simulation box, normal to the bilayer. This protein–membrane layer was solvated by 6,778 water molecules, forming a bulk liquid solution around the membrane. Cl− and K+ ions were introduced into the bulk solution to neutralize the charge of the system and establish an ionic concentration equivalent to a 150-mM solution of KCl. As the number and valency of the ions in the filter were adjusted, the net charge of the system was maintained at neutrality by adjusting the number of Cl− and K+ ions in solution.
The electrophysiological experiments of Piasta et al. (2011) used the E71A KcsA mutant to prevent their measurements of the Ba2+ blockade events from being obscured by C-type inactivation, which does not occur in the E71A mutant. This type of inactivation occurs on a much longer timescale than the length of our simulations, and no conformational change corresponding to this type of inactivation was observed in our simulations. As the E71A mutant otherwise shows the same conductivity as the wild-type channel modeled in this report, our simulations of Ba2+ blockades of the wild-type structure should generally be applicable to this mutant as well.
All simulations performed here made use of the CHARMM PARAM22 protein force field (MacKerell et al., 1998) with the backbone dihedral Cmap potential correction (MacKerell et al., 2004). The CHARMM PARAM27 force field force was used for the lipid parameters (Feller et al., 1997). Optimized Lennard-Jones parameters for Na+ and K+ were used (Noskov et al., 2004). The parameters for Ba2+ (Emin = −0.150 kcal/mol, Rmin = 1.849 Å) were determined by adjusting the Lennard-Jones radius of Ba2+ to match the experimental hydration free energy. These parameters were tested by comparing the computed radial distribution function (RDF) calculated using these parameters to a QM/MM molecular dynamics simulation (see Rowley and Roux  and Riahi et al.  for the details of this type of simulation). The RDF maxima of these two methods were identical within 0.1 Å, indicating that the Ba2+ parameters are reasonable. All other Lennard-Jones interactions were calculating using the combination rule, with the exception of the Ba2+ − O(carbonyl) interactions (Emin [Ba2+ − O] = −0.134 kcal/mol, Rmin [Ba2+ − O] = 3.36 Å), which were adjusted to reproduce the relative RIMP2/def2-TZVP hexacoordinate ion-ligand binding energies of N-methylacetamide and water. Water molecules were described using the TIP3P model (Jorgensen et al., 1983). Bonds containing hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977). The electrostatic interactions were computed with the particle mesh Ewald (PME) method, with a 72 Å × 72 Å × 81 Å grid (roughly 1 grid point per angstrom; Essmann et al., 1995). The systems were simulated with a time step of 2 fs. The temperature and pressure of the system was regulated by the CPTA method (Feller et al., 1995, 1997), where the surface area of the membrane in the xy plane was kept constant while the length of the unit cell along the z axis was allowed to vary to preserve a constant pressure.
Umbrella sampling simulations
The umbrella sampling simulations of the potential of mean force (PMF) were performed by applying planar harmonic biasing positions to the Z coordinates of the translocating ions. The details of this method are presented in a paper by Bernèche and Roux (2001). The 2D PMFs were calculated by an umbrella sampling simulation where both the Z coordinate of the Ba2+ ion (ZBa) and the Z coordinate of the lock-in ion (ZL) were simultaneous restrained. An initial grid of ion positions along the ZBa axis was defined based on the locations of the ion binding sites in the 1K4C crystal structure. For translocation of Ba2+ from site S2 to S1, this grid covered the range ZBa = 0–5.5 Å in increments of 0.5 Å. The range of ZL positions was initially assumed to extend from ZBa = 3.5 Å to 14 Å. To generate stable structures from these artificially constructed systems, energy minimizations were performed where the positions of the ions, bulk water molecules, lipids, and most of the proteins were fixed. The α-carbon of the selectivity filter residues was restrained by a 5-kcal Å−2 harmonic force and the water molecules within the channel were unrestrained. These systems were minimized by the adaptive Newton–Raphson algorithm for 200 steps. The constraints were then removed and 300-ps-long molecular dynamics simulations were performed on these systems to equilibrate the structure with the Ba2+ and lock-in ions restrained to their assigned positions. A spring constant of 10 kcal mol−1 Å−2 was used for the planar harmonic restraint on the lock-in ion, while a stronger spring constant of 50 kcal mol−1 Å−2 was used to restrain the position of the Ba2+ ion. The lock-in ion was restrained to remain inside a cylinder with a 3-Å radius around an axis running along the center of mass of the filter using a flat-bottom half harmonic restraint with a force constant of 10 kcal Å−2.
The 1D PMFs of K+ binding to the external sites of the filter when site S2 was occupied by K+ or Ba2+ presented in Fig. 2 were calculated by integrating over the regions of the corresponding 2D PMFs that correspond to Ba2+ being bound in site S2 (0 Å < ZBa < 2 Å).(1)
Hamiltonian-replica exchange MD
To improve the rate of convergence and ergodicity of our configurational sampling of the umbrella sampling windows used to calculate the PMFs, we used Hamiltonian replica exchange molecular dynamics (H-REMD) to allow exchanges between windows. In this method, attempts are made to exchange the configurations of neighboring systems on the PMF, with the acceptance of these attempts governed by the Metropolis algorithm (Jiang and Roux, 2010). This type of strategy has been demonstrated to significantly improve the convergence of distributions sampled using molecular dynamics simulations.
The 2D PMFs calculated in this study required an REMD scheme that is more elaborate than usual because umbrella-sampling windows are not distributed in a Cartesian grid, where the nearest neighbors provide a natural list of partners for exchanges. Additionally, because of the steepness of the free energy surface along the ZBa coordinate, the restraining forces on the Ba2+ ion are large (50 kcal mol−1 Å−2), leading to low exchange rates between these windows. Although reducing the spacing of the ZBa windows could improve the exchange rate, this would greatly increase the number of windows. To resolve these problems, we used a novel strategy we term “woven” H-REMD. In this method, a 1D Hamiltonian exchange sequence is constructed from the 2D grid of windows by all windows with given values of ZBa in a 1D sequence. At the maximum or minimum value of ZK for a given value of ZBa, exchanges are also attempted with a neighboring replica with an adjacent value of ZBa. This strategy allows irregularly shaped umbrella sampling grids to be sampled within the existing Hamiltonian replica exchange molecular dynamics method implemented in CHARMM. To limit the number of replicas, the windows were separated by increments of 0.5 Å, which still allowed for significant overlap and an average exchange acceptance probability of 20%. In total, the 2D PMFs for the permeation of Ba2+ from site S2 to site S1 required 162 replicas.
The unbiased 2D PMFs were computed from the umbrella sampling simulations by the weighed histogram analysis method (WHAM; Roux, 1995). The histogram bin width was 0.1 Å. The WHAM equations were iterated until convergence was achieved within a tolerance of 10−8 for all windows.
Definition of 1D translocation PMF
Let U(X) represent the potential energy as a function of all atomic coordinates X in the system, and n represent the number of ions bound to the outer sites of the pore. Let Hn(X) be an indicator function equal to 1 when the system is in state n, and zero otherwise. If we have a complete set of indicator functions, then by normalization. For the binding of a monovalent cation to the outer sites of the selectivity filter, we have only two states, with or without a bound cation. For the average of a property A(X),(2)where Pn = 〈Hn〉 is the probability the state n, and 〈A〉(n) represents the conditional average,
Let Z represent the position of the barium ion along the channel axis. For a Z-constrained average of a property A(X), the expression is a little more complicated,
Here, X′ is the integration variable and Z is the fixed constraint.
The mean force along the Z coordinate is then,
For the Ba2+ lock-in effect, we consider two states: n = 0, where there is no lock-in ion, and n = 1, where there is,where Keq(Z) is the equilibrium binding constant of the outer cation when Ba2+ is held fixed at Z. A delta function is used to restrict the configurational integrals to the holo and apo states. The equilibrium binding constant of the outer ion should be written as,
This expression is valid for any arbitrary position Z. Hence,
The free energy of substitution of K+ and Na+ (ΔΔGK→Na) was calculated using FEP molecular dynamics (FEP/MD) as the difference between the free energy to replace K+ with Na+ in a binding site and the free energy for the same transformation in bulk water (ΔΔGK→Na = ΔGK→Na,site − ΔGK→Na,bulk). The new potential energy function is defined as a linear combination of the K+ (UK) and Na+ containing (UNa) potentials according to the formula
The free energy difference of these two systems was calculated using the relation
The FEP was performed using a total of 11 values of λ in increments of 0.1 between 0 and 1. Each window was equilibrated for 200 ps before a production sampling period of 500 ps. The integration of these windows to determine the free energy difference and the sampling error was performed using the multistate Bennett acceptance ratio (MBAR) method (Shirts and Chodera, 2008). The lock-in ion was restrained to remain inside site S0 by flat-bottomed planar harmonic restraints, which have zero force when the Z coordinate of the ion is inside the site (defined by the carbonyl ligands at the top and bottom of the site) but experience a strong quadratic restraining force (krest = 100 kcal mol−1 Å−2) if the ion moves outside the site.
Calculation of the transmission coefficient
We estimated the rate of Ba2+ permeation from site S2 to site S1 when there is no lock-in ion using Grote–Hynes rate theory (Grote and Hynes, 1980). Eyring transition state theory neglects dissipative effects that diminish the reaction rate. This can be corrected for by calculating the reaction transmission coefficient, κ,
The transmission coefficient can be computed in a straightforward fashion using Grote–Hynes theory,where ω≠ is the frequency corresponding to the curvature of the PMF at the top of the barrier, i.e., , and s is the “reactive” frequency corresponding to oscillations across the barrier. This reactive frequency can be determined from the equationwhere is the Laplace transform of time-dependent friction at the top of the barrier,
can be determined by a technique proposed by Straub et al. (1988) based on analysis in terms of the generalized Langevin equation for a harmonic oscillator.(3)
Where Z is the velocity of the ion along the reaction coordinate and δZ is the displacement of the system from the transition state (δZ = Z − ZTS), and C(t) is the velocity autocorrelation function. The Laplace transform of Eq. 3 gives(4)
By substituting and , Eq. 4 can be simplified to,(5)
Using Eq. 5, can be calculated using a molecular dynamics simulation that is tightly restrained to the top of the barrier. and are calculated from ensemble averages of this trajectory. is calculated by a numerical Laplace transform of the velocity autocorrelation function from the molecular dynamics simulation.
The well frequency was determined by approximating the no lock-in PMF in Fig. 6 as a parabola with a curvature of W″ (Z╪) = 25 kcal mol−1 Å−2. To determine 〈〉TS, 〈δZ2〉TS, and the velocity autocorrelation function (C(t)), a 500-ps molecular dynamics simulation was performed where the Ba2+ ion was restrained at the top of the barrier with a 50 kcal mol−1 Å−2 harmonic restraint. These data were used to calculate numerically over a range of values of s. The reactive frequency, s, was determined through a numerical solution of Eq. 4 using the value of ω╪ from the PMF and these values of , yielding a reactive frequency of s = 3 ps−1 and a transmission coefficient of κ = 0.25.
RESULTS AND DISCUSSION
PMF of lock-in ion binding
Piasta et al. (2011) used the available structural and electrophysiological data for KcsA to assign the Ba2+ binding and lock-in sites. Analysis of the K+-free block time distribution was consistent with two distinct Ba2+ binding sites within the filter. Because the electronic density associated with Ba2+ was only observed in sites S2 and S4 in the crystal structure of KcsA soaked with Ba2+ (Lockless et al., 2007), it was assumed that sites S2 and S4 are the sites that are predominantly occupied by Ba2+ during the blocking events. Based on the analysis of electrophysiological data, it was inferred that the long blocks occur predominantly when the Ba2+ is bound to site S2, implying that the external cation must be binding to either site S0 or site S1.
Piasta et al. (2011) argued that the high selectivity displayed in the experiments is inconsistent with S0 being the lock-in site, and postulated that site S1 must be the actual lock-in site. One may note that this leads to a multi-ion binding arrangement that differs from the generally accepted view, in that two ions are simultaneously bound to adjacent sites (S2 and S1). X-ray crystallography of ion-bound channels (Zhou and MacKinnon, 2003), streaming potential measurements (Alcayaga et al., 1989), and molecular dynamics simulations (Bernèche and Roux, 2000) have all indicated that alkali ions preferentially occupy alternating binding sites within the filter, with water molecules interspersed in between. This arrangement avoids the strong electrostatic repulsion that would occur when ions occupy neighboring sites, an effect that intuitively should be even larger for the interaction between the divalent Ba2+ ion and an alkali lock-in ion.
To evaluate the proposed scenario within the framework of molecular dynamics (MD) simulations, we calculated the PMF of the lock-in K+ ion moving along the Z axis of the filter from site S1 to the external solution when there is a Ba2+ ion occupying site S2 (Fig. 2). For comparison, we also show the PMF for this ion when site S2 is occupied by a K+ ion. The PMF between ZK = 3 Å and ZK = 5 Å corresponds to the K+ ion occupying site S1. This is a highly unstable arrangement when site S2 is occupied by Ba2+, with these configurations lying >10 kcal/mol higher on the free energy surface. This configuration is far less stable in comparison to the configuration where sites S2 and S4 are occupied by K+ ions, in keeping with the higher electrostatic repulsion the lock-in ion experiences with a divalent Ba2+ ion in comparison to a K+ ion. This effect diminishes when the lock-in K+ ion is in a more external binding site; the binding energy of site S0 is comparable for the Ba2+ and the K+ occupied filters, although the lock-in K+ ion tends to bind more outwardly in site S0 when Ba2+ is present. These PMFs show that K+ is unlikely to bind in site S1 when site S2 is occupied by Ba2+. A scenario where site S2 is the external Ba2+ binding site and S1 is the K+ binding lock-in site is ruled out based on these considerations.
PMF calculations of barium permeation
We considered the alternative scenario where K+ binding to site S0 is responsible for blocking the forward translocation of Ba2+ from site S2 to S1. Previous computational studies found that site S0 is not selective for K+ because it is wider than the lower sites, and ions bound in it are partially exposed to the external solution (Noskov et al., 2004; Noskov and Roux, 2006); however, a computational study by Kim and Allen (2011) found that site S0 is selective for K+ when the filter is occupied by two Ba2+ ions. To comprehensively evaluate if K+ bound in site S0 can block the permeation of Ba2+, we computed the 2D PMF of Ba2+ translocating from site S2 to S1 when there is a K+ ion external to it. To evaluate the ion selectivity of this lock-in effect, we also computed the PMF when the external ion is Na+. In each case, an unrestrained K+ ion is present in the bottom of the filter, occupying site S4 or S3, depending on the configuration of the other ions.
The 2D PMFs for the translocation of Ba2+ from S2 to S1 with an external lock-in ion are presented in Fig. 3. In these 2D free energy maps, the x axis corresponds to the position of the Ba2+, moving from site S2 (0 Å < ZBa < 2 Å) to site S1 (3 Å < ZBa < 5 Å), while the y axis corresponds to the position of the lock-in ion, K+ (left) or Na+ (right). All ion positions are defined relative to the center of mass of the filter. The PMF for the translocation of Ba2+ from S2 to S1 shows three distinct minima corresponding to metastable configurations: Ba2+ bound in S2 with the lock-in ion in S0 (Fig. 3, bottom left), the Ba2+ bound in S2 with the lock-in ion is the external solution (Fig. 3, top left), and Ba2+ bound in S1 with the lock-in ion in the external solution (Fig. 3, top right).
The physical basis for the lock-in effect is immediately apparent from these free energy surfaces; the barriers are extremely high for direct translocation when the lock-in ion is present in S0 (ZL > 0 Å). This reflects that the translocation of Ba2+ to site S1 while K+ occupies S0 leads to an arrangement where the ions occupy neighboring sites and experience a strongly repulsive electrostatic interaction. The lowest free energy path for Ba2+ translocation corresponds to the complete exit of the lock-in ion into the external solution followed by the translocation of Ba2+ from S2 to S1. The lowest barrier for this transition relative to the minimum is 17 kcal/mol, which is consistent with a slow translocation rate of Ba2+ in comparison to alkali ions. After crossing this barrier, Ba2+ is bound in site S1 and the lower unrestrained ion has progressed to site S3. The free energy of the state where Ba2+ is bound in S2 is within 1 kcal/mol of the state where Ba2+ is bound in S1. This indicates that sites S1 and S2 have similar affinities for Ba2+ in the multi-ion binding scenario. However, once the Ba2+ is bound in S1, external ions can no longer impede its progress toward the exit of the pore, as there are no binding sites beyond S0. For this reason, within this model, the lock-in effect reflects the critical attempt of Ba2+ to translocate from S2 to S1.
The global energy minimum of these surfaces (near ZBa = 3.5 Å and ZL > 12 Å) corresponds to Ba2+ occupying site S1 and the lock-in ion in the external solution. The energy minimum of the state where Ba2+ occupies site S2 and K+ occupies site S0 is 7 kcal/mol higher in energy. This may be an underestimate of the binding affinity of site S0, as the nonpolarizable force field models used in these studies can underestimate the interaction between ions and carbonyls (Yu et al., 2010).
The PMFs of the Na+ and K+ lock-in ions are very similar when the ions are outside of site S0, which is consistent with a scenario in which the ion is in the external solution and therefore effectively uncoupled from the filter. The interesting distinctions occur when the lock-in ion is in S0 (ZL < 9 Å) and Ba2+ is in S2 (ZBa < 2 Å). The minimum on K+ PMF occurs when the ion occupies the middle of S0, near ZL = 7.5 Å. The minimum of the PMF for Na+ occurs lower, at ZL = 7 Å, where it can interact closely to the four backbone carbonyls. This portion of the PMF is ∼2 kcal/mol higher in energy for Na+ than it is for K+, which is consistent with site S0 acting as a K+ selective site.
A recent computational study by Kim and Allen (2011) also calculated the PMF of the K+ and Na+ ions in positions between site S0 and the external solution while the filter was occupied by Ba2+, and proposed that the lock-in phenomenon results from K+ binding in site S0 while Ba2+ binds in site S2. Our results are generally consistent with theirs; however, it should be noted that their simulations differs from ours in that two Ba2+ ions were included simultaneously in the S4 and S2 sites in their simulations. Although Ba2+ ion densities were observed in both these sites in the structure obtained after growing the crystals in 5 mM BaCl2 (Lockless et al., 2007), Piasta et al. (2011) determined that the channel is occupied by only one Ba2+ ion under the conditions of their blockade experiments (from the K+-free block time distribution). Based on these observations, the models used in our calculations were constructed in such a way that there is only one Ba2+ ion in the filter.
FEP studies of site selectivity
A comparison of the computed PMFs with the K+ and Na+ lock-in ions indicates that site S0 is moderately selective for K+ over Na+ when Ba2+ is bound in site S2. This is somewhat surprising, as previous computational studies found that site S0 is weakly selective for Na+ when the lower sites of the filter were occupied by K+ (Noskov and Roux, 2006; Egwolf and Roux, 2010). To resolve these differences, we quantified the thermodynamic selectivity of site S0 when K+, or Ba2+, occupies S2. This was achieved by calculating ΔΔGK→Na, which is the free energy of substituting a bound K+ ion for a Na+ ion relative to this substitution in bulk water.
The free energy of selectivity of a site, ΔΔGK→Na, is difficult to determine accurately using PMFs, as the resolution of the free energy surfaces is limited and sampling error can cause sizable deviations. Furthermore, this method relies on a stable baseline corresponding to the ion in the bulk, which requires extending the PMF for the ion reaching a large distance away from the pore. Computationally, it is more advantageous and straightforward to calculate the free energy of this exchange using alchemical FEP. This method computes the free energy of converting from one potential energy function to another. In this instance, we are computing the free energy difference of replacing the K+ ion in site S0 with a Na+ ion. We performed two FEP simulations: one where site S2 is occupied by a Ba2+ ion and one where it is occupied by a K+ ion. These free energies of selectivity differ from those reported in the previous studies by Noskov and Roux (2006) because we use flat bottom restraints and the simulation temperature used here is 25°C, which is the temperature used in the blockade experiments of Piasta et al. (2011), as opposed to Noskov and Roux (2006), where the physiological temperature of 37°C was used. The ion filter configurations and the corresponding free energy of selectivity are illustrated in Fig. 4. The details of these calculations are presented in the Materials and methods section.
Consistent with previous studies, site S0 is essentially nonselective when S2 and S4 are occupied by K+, with a ΔΔGK→Na = −0.5 kcal/mol. Remarkably, this value changes significantly when a Ba2+ ion occupies site S2, leading site S0 to become modestly selective for K+, with a ΔΔGK→Na = +1.8 kcal/mol. This is consistent with the 2D PMFs presented in the “PMF of lock-in ion binding” section, which are ∼2 kcal/mol higher for Na+ in the areas corresponding to the lock-in ion occupying site S0. Therefore, the presence of Ba2+ is responsible for a shift of ∼2.3 kcal/mol.
To understand the origin of this difference, we calculated the axial distribution function of the lock-in ions when bound in site S0 (Fig. 5). When K+ occupies site S2, Na+ occupies a broad range of positions within the site, but tends to bind at the inner edge of S0, where it can coordinate directly to the four backbone carbonyls, a feature noted previously (Shrivastava et al., 2002; Kim and Allen, 2011), while K+ tends to bind in a more outer fashion in the site (further away toward the extracellular side), coordinating with both the Tyr78 carbonyl ligands and water molecules from the external solution. When Ba2+ occupies S2, both ions are pushed outward, with most probable positions near Z = 7.5 Å. But the impact on the coordination environment of K+ is minimal. In contrast, Na+ is shifted to a significantly more outer position; its coordination environment becomes more similar to that of K+, which causes this site to become K+ selective when Ba2+ occupies site S2. The trend from the free energy simulations is in qualitative agreement with the experimental estimates of Piasta et al. (2011). However, there is a considerable quantitative discrepancy, as the calculated ΔΔGK→Na of +1.8 kcal/mol is considerably smaller than the experimental estimate of >7 kcal/mol.
PMF interpretation of barium blockades
The 2D PMFs display features that are qualitatively consistent with the concept of a Ba2+ blockade that is affected by the presence of an external lock-in ion, although this description does not directly demonstrate how the rate of Ba2+ exit into the external solution corresponds to these PMFs. It is useful to clarify this matter further by calculating the rate of permeation of the Ba2+ ion from site S2 to S1, which can be described using transition state theory (Glasstone et al., 1941; Chandler, 1978),(6)where W is the PMF corresponding to the translocation of the Ba2+ outwards from the filter along reaction coordinate Z, ΔW╪ is the height of the PMF relative to the reactant well, Zm is the location of the minimum corresponding to Ba2+ occupying site S2, Z╪ is the location of the maximum, and κ is the transmission coefficient (0 < κ < 1). The calculation of these terms is described in the Materials and methods section.
Although Eq. 1 is a familiar expression, the effect of the lock-in ion is implicitly included in the activation free energy ΔW╪. This is where the difficulty lies because constructing the generalized 1D PMF for the permeation of Ba2+ is more complicated, as the PMF depends on the concentration of the external cation. As the 2D PMFs illustrated, the presence of the lock-in ion can considerably affect the 1D PMFs corresponding to Ba2+ permeation. The probability of an ion being bound to the filter depends on the external concentration of the lock-in ion, so application of transition state theory requires that we define a 1D PMF of the ZBa coordinate that varies as a function of concentration of the lock-in ion in the external solution. To this end, we express the mean force as a linear combination of the mean force when there is no external lock-in ion, Wn=0(ZBa) and the PMF when there is an external lock-in ion Wn=1(ZBa). These mean forces are weighted by the probability (P) of the lock-in ion being present when the Ba2+ is present at point ZBa.
The probabilities can be expressed in terms of the concentration of the lock-in ion in the external solution ([C]) and the equilibrium constant of binding to the lock-in site when Ba2+ is at ZBa (Keq(ZBa)),
The derivation of this equation is presented in the Materials and methods section.
In principle, a PMF along the Z axis could be determined by integrating Eq. 2, although for our purposes, we only need to consider the two extreme cases: (1) the [C] = 0 case, when there are no blocking ions in the external solution and the only contribution is from , and (2) when [C] is at a saturating concentrations and will dominate. By calculating the PMF for these two limiting cases, we can describe the lock-in effect.
We calculated the 1D PMFs of Ba2+ permeation subject with lock-in ions Na+ or K+ external to it from the 2D PMF of Ba2+. The 1D PMF of Ba2+ permeation when the lock-in ion is in site S0, Wn=1(ZBa), can be estimated by integrating over the regions of the 2D PMF that correspond to the lock-in ion being present in site S0 (ZL = [5 Å, 9 Å]),
The PMF when there is no lock-in ion present can be estimated by averaging over the regions of the 2D PMF where the lock-in ion is too distant to affect permeation (ZL = [9 Å, 14 Å]),
We averaged the Na+ and K+ PMF to remove spurious differences due to sampling error. The 1D PMF for this unblocked permeation is presented in Fig. 6. When site S0 is not occupied by an ion, permeation of Ba2+ occurs with a barrier of roughly 15 kcal/mol. The Eyring reaction rate calculated using Eq. 6 for Ba2+ permeation when there is no lock-in ion present is 208 s−1, and the Grote–Hynes transmission coefficient is κ = 0.25, yielding a rate constant of ∼52 s−1, which is in reasonable agreement with the experimentally determined rate constant of 204 s−1. This is consistent with a slow but nonzero rate of permeation of Ba2+ in the absence of a lock-in ion.
The second limiting case, where ion concentration in the external solution is sufficiently high to saturate site S0, is plotted in Fig. 6 in purple for Na+ and green for K+. There are minor differences between the Na+ and K+ curves; however, the more significant feature is that both PMFs are effectively infinite except when Ba2+ is in its initial state in site S2. This reflects the high repulsion between the two ions when they are placed in neighboring sites. In this regime, the rate of outward permeation will be zero for both types of lock-in ions. This corresponds to the kinetic model of Piasta et al. (2011), which determined that the rate of outward permeation of Ba2+ is effectively zero when the external concentration of K+ was >1 mM.
Within this framework, the strong lock-in effect observed in the experiments of Miller and coworkers (Neyton and Miller, 1988; Piasta et al., 2011) would occur if the equilibrium constant for an ion binding to the lock-in site were large. A K+ selective lock-in effect would occur if the equilibrium constant were larger for K+ than Na+. Our computations are not quantitatively consistent with the high selectivity and micromolar affinity for K+ that was experimentally observed by Piasta et al. (2011). The absolute binding affinity of K+ to site S0 is is too low according to the calculated PMF (Fig. 2) and and the site is only weakly selective for K+ over Na+ when site S2 is occupied by Ba2+ according to FEP (Fig. 5).
Possible explanations for the lack of quantitative agreement
Although the present simulations capture some important features of a selective lock-in effect of a Ba2+ blockade, they are not quantitatively consistent with the experiment. Such inconsistencies can be attributed to the limitations of various aspects of the computational treatment. For instance, as we used a nonpolarizable force field, the induced electronic polarization of the environment by the ions is neglected in these models. This is of particular concern when considering the effect of the medium screening the interaction between a divalent Ba2+ ion and the lock-in ion. The development of a more accurate force field accounting explicitly for the effect of induced polarization will, hopefully, allow a more satisfactory representation of ion permeation (Lamoureux et al., 2003; Lopes et al., 2007; Yu et al., 2010; Rowley and Roux, 2012). Furthermore, we used the high-resolution crystallographic structure of the KcsA channel (PDB accession no. 1K4C), which corresponds to a functional state with a closed intracellular gate. Although the processes of interest are taking place in the selectivity filter away from the intracellular gate, a high-resolution x-ray structure of the channel in the open-conductive state might provide a more realistic system for the simulations. Lastly, selectivity is governed by small free energy differences, and achieving statistically converged results is challenging, though feasible through enhanced sampling methods. Although we can attribute some of the discrepancies between our simulations to these technical limitations, current explanations for the lock-in effect might also need to be examined. In particular, the calculated PMFs show that stability of Ba2+ binding in site S1 is comparable to its binding in site S2.
The possibility that Ba2+ might bind to sites in the selectivity filter other than sites S2 and S4 will require further consideration. Electronic density from the Ba2+ ions was observed only for sites S2 and S4 in the crystallographic structure, whereas no density was observed for sites S1 and S3 (Lockless et al., 2007). The resolution of the Ba2+-occupied x-ray structure is too low to unambiguously determine the occupancy of the binding sites, so there are several plausible binding configurations that ought to be considered in analyzing the functional data. One interpretation, adopted by Piasta et al. (2011), is that Ba2+ does not bind favorably in sites S1 and S3. Assuming that ∼10% occupancy is a reasonable minimum threshold for detection of ion binding, this would imply that the energy of Ba2+ in sites S1 and S3 is less favorable by at least 1.4 kcal/mol relative to sites S2 or S4. Thus, a very small energy difference could explain the binding pattern observed in the x-ray structure. Alternatively, the observed electronic density might reflect a structure where two Ba2+ ions are bound simultaneously in sites S2 and S4. The x-ray structure was obtained from crystals that were grown from a solution with 5 mM BaCl2, which is considerably larger than the micromolar Ba2+ concentrations used in the single-Ba2+ occupancy blockade experiments (Piasta et al., 2011). The implication is that single occupancy of site S3 may be energetically accessible, though not observed in the x-ray structure due to the high electrostatic penalty associated with divalent ions occupying adjacent sites. Greater clarity about the occupancy of the Ba2+ in the binding sites is essential to definitively interpreting these experiments.
The simulations in this study were undertaken to help interpret the Ba2+ block experiments of the KcsA channel performed by Piasta et al. (2011). These calculations were designed to examine the translocation of the Ba2+ ion from site S2 to site S1 and how this specific process is affected by the binding of extracellular cations to the filter. This choice was motivated by two observations. First, the Ba2+ ions in the Ba2+-soaked crystal structure of KcsA bind either to site S2 or site S4 of the selectivity filter. Second, the electrophysiological data suggested that the dominant situation during the prolonged blocks occurs when Ba2+ is bound to the outermost of these two sites. Under these conditions, the external lock-in site must either be site S1 or site S0, two possibilities that were investigated computationally.
Using umbrella sampling PMF calculations, we determined the affinity of K+ to sites S0 and S1 when site S2 is occupied by either Ba2+ or K+. The computations indicate that binding of a cation (K+ or Na+) in site S1 while a Ba2+ is bound in site S2 is energetically prohibitive. Although site S1 has a weak affinity for K+ when S2 is occupied by K+, it is highly unfavorable for binding K+ when S2 is occupied by Ba2+. This suggests that the proposal of Piasta et al. (2011) that the lock-in effect is due to a block of the translocation of Ba2+ from site S2 to site S1 caused by the binding of a K+ ion in the adjacent site S1 is unlikely.
If S1 is not the lock-in site, the only remaining external site that could serve in this capacity is S0. This conclusion is somewhat unexpected and counterintuitive. Previous computational studies have indicated that S0 is not selective for K+. In fact, this observation partly motivated Piasta et al. (2011) to rule out the possibility that S0 may be able to act as a lock-in site for the selective binding of external cations. However, the PMF calculations presented in this paper indicate that S0 becomes considerably more selective when S2 is occupied by Ba2+ than had been observed in previous studies when the filter was occupied by K+. The trend was confirmed by FEP calculations indicating that site S0 sites became significantly more selective for K+ when Ba2+ is bound inside the filter, although it is worth noting that the computed selectivities are smaller than the values determined by the experiments of Piasta et al. (2011). We attribute the increase in selectivity to the greater electrostatic repulsion between the divalent Ba2+ ion and the lock-in ion. This repulsion drives the lock-in ion bound to site S0 in a more outward position, which disproportionately affects Na+, as it tends to bind favorably in a more inward position in the site. This suggests that the selectivity derived from Ba2+ blockade experiments differs from the selectivity that occurs under physiological conditions.
To clarify the mechanism of block by the external lock-in site, we calculated the rate of Ba2+ permeation using Grote–Hynes theory and our calculated 1D PMF of Ba2+ translocation from S2 to S1. The calculated rate constant for this process is ∼52 s−1, which is in reasonable agreement with the experimental rate of 204 s−1. A statistical mechanical formulation of the rate was elaborated to incorporate the concentration dependence of the external lock-in ion. When there is a saturating concentration of the lock-in ion in the external solution, the PMF for permeation has an extremely high barrier. This reflects a scenario where Ba2+ permeates from S2 to S1 when there is an alkali ion in S0. The weighting of these PMFs depends on the equilibrium constant for the binding of the lock-in ion from the external solution, where a high binding affinity will increase the weighting of the “locked-in” PMF. Nevertheless, although the calculated Ba2+ permeation PMFs and estimated rate constant were consistent with this kind of lock-in scenario, the affinity of K+ to site S0 and the small selectivity of K+ over Na+ in this site are not in quantitative agreement with the experiment.
The lack of quantitative agreement between the calculations and the experimental data is not surprising and confirms that computational models are imperfect and need to be improved. Various limitations and deficiencies, such as sampling error, inaccuracy of nonpolarizable force fields, and the appropriateness of the crystallographic structures used in the simulations and how they correspond to the functional data should also be considered. Although the present study was limited to an examination of the translocation of the Ba2+ ion from site S2 to site S1, it is possible that the observed lock-in effects actually involve translocation steps with Ba2+ bound to sites that are more inward (e.g., S3 or S4). Ultimately, this will require a complete characterization of each translocation step involved when a Ba2+ ion enters the filter from the intracellular side and proceeds all the way across the pore to exit on the extracellular side. In our view, such an ambitious undertaking would be premature at this point given the current limitations of the force field and computational model. Gradual progress is being made in these areas, which provides the prospect that we will eventually be able to achieve quantitative agreement between this experimental electrophysiological data and molecular dynamics simulations.
We thank Dr. Balasundaresan Dhakshnamoorthy for assistance in interpreting the crystallographic structural data and Jessica Besaw for a thorough proofreading of the manuscript.
This work was supported by the National Institutes of Health (NIH) through grant NIH/NIGMS R01-GM062342. C.N. Rowley acknowledges the Natural Sciences and Engineering Research Council of Canada for a postdoctoral fellowship and a discovery grant. The simulations reported in this paper were performed using computational resources at the Laboratory Computing Resource Center at the Argonne National Laboratory, the University of Chicago Argonne Computation Institute, the ACEnet and SciNet facilities of the Compute Canada consortium (RAPI: djk-615-ab), and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. OCI-1053575.
Christopher Miller served as editor.
- Abbreviations used in this paper:
- free energy perturbation
- potential of mean force
- Submitted: 21 June 2013
- Accepted: 16 August 2013
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/).