A computational study of barium blockades in the KcsA potassium channel based on multi-ion potential of mean force calculations and free energy perturbation

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 S₂, 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ ions was introduced in the “internal” solution, the single channel current underwent frequent interruptions due to Ba²⁺ binding to the selectivity filter. The block is relieved when the Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ blockade experiments are highly complementary to detailed free energy computations using molecular dynamics simulations. As the permeation of Ba²⁺ 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 Ba²⁺ between sites. The Ba²⁺ 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 Ba²⁺ blockade experiments of Piasta et al. (2011) in an attempt to reconcile the experimental observations with the results of simulation.

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 Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ (E_min = −0.150 kcal/mol, R_min = 1.849 Å) were determined by adjusting the Lennard-Jones radius of Ba²⁺ 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 [2012] and Riahi et al. [2013] for the details of this type of simulation). The RDF maxima of these two methods were identical within 0.1 Å, indicating that the Ba²⁺ parameters are reasonable. All other Lennard-Jones interactions were calculating using the combination rule, with the exception of the Ba²⁺ − O(carbonyl) interactions (E_min [Ba²⁺ − O] = −0.134 kcal/mol, R_min [Ba²⁺ − 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.

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 Ba²⁺ ion (Z_Ba) and the Z coordinate of the lock-in ion (Z_L) were simultaneous restrained. An initial grid of ion positions along the Z_Ba axis was defined based on the locations of the ion binding sites in the 1K4C crystal structure. For translocation of Ba²⁺ from site S₂ to S₁, this grid covered the range Z_Ba = 0–5.5 Å in increments of 0.5 Å. The range of Z_L positions was initially assumed to extend from Z_Ba = 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 Å⁻² 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 Ba²⁺ and lock-in ions restrained to their assigned positions. A spring constant of 10 kcal mol⁻¹ Å⁻² was used for the planar harmonic restraint on the lock-in ion, while a stronger spring constant of 50 kcal mol⁻¹ Å⁻² was used to restrain the position of the Ba²⁺ 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 Å⁻².

The 1D PMFs of K⁺ binding to the external sites of the filter when site S₂ was occupied by K⁺ or Ba²⁺ presented in Fig. 2 were calculated by integrating over the regions of the corresponding 2D PMFs that correspond to Ba²⁺ being bound in site S₂ (0 Å < Z_Ba < 2 Å).

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 Z_Ba coordinate, the restraining forces on the Ba²⁺ ion are large (50 kcal mol⁻¹ Å⁻²), leading to low exchange rates between these windows. Although reducing the spacing of the Z_Ba 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 Z_Ba in a 1D sequence. At the maximum or minimum value of Z_K for a given value of Z_Ba, exchanges are also attempted with a neighboring replica with an adjacent value of Z_Ba. 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 Ba²⁺ from site S₂ to site S₁ 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⁻⁸ for all windows.

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 H_n(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 $∑nHn=1$ 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),

where P_n = 〈H_n〉 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 Ba²⁺ lock-in effect, we consider two states: n = 0, where there is no lock-in ion, and n = 1, where there is,

where K_eq(Z) is the equilibrium binding constant of the outer cation when Ba²⁺ 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⁺ (ΔΔG_K→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 (ΔΔG_K→Na = ΔG_K→Na,site − ΔG_K→Na,bulk). The new potential energy function is defined as a linear combination of the K⁺ (U_K) and Na⁺ containing (U_Na) 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 S₀ 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 (k_rest = 100 kcal mol⁻¹ Å⁻²) if the ion moves outside the site.

We estimated the rate of Ba²⁺ permeation from site S₂ to site S₁ 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., $ω╪2=‖W″(Z╪)‖/m$⁠, and s is the “reactive” frequency corresponding to oscillations across the barrier. This reactive frequency can be determined from the equation

where $M^(s)$ is the Laplace transform of time-dependent friction at the top of the barrier,

$M^(s)$ 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.

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 − Z_TS), and C(t) is the velocity autocorrelation function. The Laplace transform of Eq. 3 gives

By substituting $C(0)=〈Z.2〉TS$ and $m=kBT/〈Z.2〉TS$⁠, Eq. 4 can be simplified to,

Using Eq. 5, $M^(s)/m$ can be calculated using a molecular dynamics simulation that is tightly restrained to the top of the barrier. $〈Z.2〉TS$ and $〈δZ2〉TS$ are calculated from ensemble averages of this trajectory. $C^(s)$ 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⁻¹ Å⁻². To determine 〈$Z.2$〉_TS, 〈δZ²〉_TS, and the velocity autocorrelation function (C(t)), a 500-ps molecular dynamics simulation was performed where the Ba²⁺ ion was restrained at the top of the barrier with a 50 kcal mol⁻¹ Å⁻² harmonic restraint. These data were used to calculate $M^(s)$ 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 $M^(s)$⁠, yielding a reactive frequency of s = 3 ps⁻¹ and a transmission coefficient of κ = 0.25.

Piasta et al. (2011) used the available structural and electrophysiological data for KcsA to assign the Ba²⁺ binding and lock-in sites. Analysis of the K⁺-free block time distribution was consistent with two distinct Ba²⁺ binding sites within the filter. Because the electronic density associated with Ba²⁺ was only observed in sites S₂ and S₄ in the crystal structure of KcsA soaked with Ba²⁺ (Lockless et al., 2007), it was assumed that sites S₂ and S₄ are the sites that are predominantly occupied by Ba²⁺ during the blocking events. Based on the analysis of electrophysiological data, it was inferred that the long blocks occur predominantly when the Ba²⁺ is bound to site S₂, implying that the external cation must be binding to either site S₀ or site S₁.

Piasta et al. (2011) argued that the high selectivity displayed in the experiments is inconsistent with S₀ being the lock-in site, and postulated that site S₁ 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 (S₂ and S₁). 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 Ba²⁺ 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 S₁ to the external solution when there is a Ba²⁺ ion occupying site S₂ (Fig. 2). For comparison, we also show the PMF for this ion when site S₂ is occupied by a K⁺ ion. The PMF between Z_K = 3 Å and Z_K = 5 Å corresponds to the K⁺ ion occupying site S₁. This is a highly unstable arrangement when site S₂ is occupied by Ba²⁺, 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 S₂ and S₄ are occupied by K⁺ ions, in keeping with the higher electrostatic repulsion the lock-in ion experiences with a divalent Ba²⁺ 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 S₀ is comparable for the Ba²⁺ and the K⁺ occupied filters, although the lock-in K⁺ ion tends to bind more outwardly in site S₀ when Ba²⁺ is present. These PMFs show that K⁺ is unlikely to bind in site S₁ when site S₂ is occupied by Ba²⁺. A scenario where site S₂ is the external Ba²⁺ binding site and S₁ is the K⁺ binding lock-in site is ruled out based on these considerations.

We considered the alternative scenario where K⁺ binding to site S₀ is responsible for blocking the forward translocation of Ba²⁺ from site S₂ to S₁. Previous computational studies found that site S₀ 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 S₀ is selective for K⁺ when the filter is occupied by two Ba²⁺ ions. To comprehensively evaluate if K⁺ bound in site S₀ can block the permeation of Ba²⁺, we computed the 2D PMF of Ba²⁺ translocating from site S₂ to S₁ 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 S₄ or S₃, depending on the configuration of the other ions.

The 2D PMFs for the translocation of Ba²⁺ from S₂ to S₁ 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 Ba²⁺, moving from site S₂ (0 Å < Z_Ba < 2 Å) to site S₁ (3 Å < Z_Ba < 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 Ba²⁺ from S₂ to S₁ shows three distinct minima corresponding to metastable configurations: Ba²⁺ bound in S₂ with the lock-in ion in S₀ (Fig. 3, bottom left), the Ba²⁺ bound in S₂ with the lock-in ion is the external solution (Fig. 3, top left), and Ba²⁺ bound in S₁ 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 S₀ (Z_L > 0 Å). This reflects that the translocation of Ba²⁺ to site S₁ while K⁺ occupies S₀ leads to an arrangement where the ions occupy neighboring sites and experience a strongly repulsive electrostatic interaction. The lowest free energy path for Ba²⁺ translocation corresponds to the complete exit of the lock-in ion into the external solution followed by the translocation of Ba²⁺ from S₂ to S₁. The lowest barrier for this transition relative to the minimum is 17 kcal/mol, which is consistent with a slow translocation rate of Ba²⁺ in comparison to alkali ions. After crossing this barrier, Ba²⁺ is bound in site S₁ and the lower unrestrained ion has progressed to site S₃. The free energy of the state where Ba²⁺ is bound in S₂ is within 1 kcal/mol of the state where Ba²⁺ is bound in S₁. This indicates that sites S₁ and S₂ have similar affinities for Ba²⁺ in the multi-ion binding scenario. However, once the Ba²⁺ is bound in S₁, external ions can no longer impede its progress toward the exit of the pore, as there are no binding sites beyond S₀. For this reason, within this model, the lock-in effect reflects the critical attempt of Ba²⁺ to translocate from S₂ to S₁.

The global energy minimum of these surfaces (near Z_Ba = 3.5 Å and Z_L > 12 Å) corresponds to Ba²⁺ occupying site S₁ and the lock-in ion in the external solution. The energy minimum of the state where Ba²⁺ occupies site S₂ and K⁺ occupies site S₀ is 7 kcal/mol higher in energy. This may be an underestimate of the binding affinity of site S₀, 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 S₀, 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 S₀ (Z_L < 9 Å) and Ba²⁺ is in S₂ (Z_Ba < 2 Å). The minimum on K⁺ PMF occurs when the ion occupies the middle of S₀, near Z_L = 7.5 Å. The minimum of the PMF for Na⁺ occurs lower, at Z_L = 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 S₀ 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 S₀ and the external solution while the filter was occupied by Ba²⁺, and proposed that the lock-in phenomenon results from K⁺ binding in site S₀ while Ba²⁺ binds in site S₂. Our results are generally consistent with theirs; however, it should be noted that their simulations differs from ours in that two Ba²⁺ ions were included simultaneously in the S₄ and S₂ sites in their simulations. Although Ba²⁺ ion densities were observed in both these sites in the structure obtained after growing the crystals in 5 mM BaCl₂ (Lockless et al., 2007), Piasta et al. (2011) determined that the channel is occupied by only one Ba²⁺ 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 Ba²⁺ ion in the filter.

A comparison of the computed PMFs with the K⁺ and Na⁺ lock-in ions indicates that site S₀ is moderately selective for K⁺ over Na⁺ when Ba²⁺ is bound in site S₂. This is somewhat surprising, as previous computational studies found that site S₀ 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 S₀ when K⁺, or Ba²⁺, occupies S₂. This was achieved by calculating ΔΔG_K→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, ΔΔG_K→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 S₀ with a Na⁺ ion. We performed two FEP simulations: one where site S₂ is occupied by a Ba²⁺ 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 S₀ is essentially nonselective when S₂ and S₄ are occupied by K⁺, with a ΔΔG_K→Na = −0.5 kcal/mol. Remarkably, this value changes significantly when a Ba²⁺ ion occupies site S₂, leading site S₀ to become modestly selective for K⁺, with a ΔΔG_K→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 S₀. Therefore, the presence of Ba²⁺ 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 S₀ (Fig. 5). When K⁺ occupies site S₂, Na⁺ occupies a broad range of positions within the site, but tends to bind at the inner edge of S₀, 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 Ba²⁺ occupies S₂, 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 Ba²⁺ occupies site S₂. 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 ΔΔG_K→Na of +1.8 kcal/mol is considerably smaller than the experimental estimate of >7 kcal/mol.

The 2D PMFs display features that are qualitatively consistent with the concept of a Ba²⁺ blockade that is affected by the presence of an external lock-in ion, although this description does not directly demonstrate how the rate of Ba²⁺ 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 Ba²⁺ ion from site S₂ to S₁, which can be described using transition state theory (Glasstone et al., 1941; Chandler, 1978),

where W is the PMF corresponding to the translocation of the Ba²⁺ outwards from the filter along reaction coordinate Z, ΔW^╪ is the height of the PMF relative to the reactant well, Z_m is the location of the minimum corresponding to Ba²⁺ occupying site S₂, 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 Ba²⁺ 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 Ba²⁺ 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 Z_Ba 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, W_n=0(Z_Ba) and the PMF when there is an external lock-in ion W_n=1(Z_Ba). These mean forces are weighted by the probability (P) of the lock-in ion being present when the Ba²⁺ is present at point Z_Ba.

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 Ba²⁺ is at Z_Ba (K_eq(Z_Ba)),

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 $〈F〉(ZBa,n=0)$⁠, and (2) when [C] is at a saturating concentrations and $〈F〉(ZBa,n=1)$ will dominate. By calculating the PMF for these two limiting cases, we can describe the lock-in effect.

We calculated the 1D PMFs of Ba²⁺ permeation subject with lock-in ions Na⁺ or K⁺ external to it from the 2D PMF of Ba²⁺. The 1D PMF of Ba²⁺ permeation when the lock-in ion is in site S₀, W_n=1(Z_Ba), can be estimated by integrating over the regions of the 2D PMF that correspond to the lock-in ion being present in site S₀ (Z_L = [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 (Z_L = [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 S₀ is not occupied by an ion, permeation of Ba²⁺ occurs with a barrier of roughly 15 kcal/mol. The Eyring reaction rate calculated using Eq. 6 for Ba²⁺ permeation when there is no lock-in ion present is 208 s⁻¹, and the Grote–Hynes transmission coefficient is κ = 0.25, yielding a rate constant of ∼52 s⁻¹, which is in reasonable agreement with the experimentally determined rate constant of 204 s⁻¹. This is consistent with a slow but nonzero rate of permeation of Ba²⁺ 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 S₀, 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 Ba²⁺ is in its initial state in site S₂. 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 Ba²⁺ 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 S₀ 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 S₂ is occupied by Ba²⁺ according to FEP (Fig. 5).

Although the present simulations capture some important features of a selective lock-in effect of a Ba²⁺ 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 Ba²⁺ 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 Ba²⁺ binding in site S₁ is comparable to its binding in site S₂.

The possibility that Ba²⁺ might bind to sites in the selectivity filter other than sites S₂ and S₄ will require further consideration. Electronic density from the Ba²⁺ ions was observed only for sites S₂ and S₄ in the crystallographic structure, whereas no density was observed for sites S₁ and S₃ (Lockless et al., 2007). The resolution of the Ba²⁺-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 Ba²⁺ does not bind favorably in sites S₁ and S₃. Assuming that ∼10% occupancy is a reasonable minimum threshold for detection of ion binding, this would imply that the energy of Ba²⁺ in sites S₁ and S₃ is less favorable by at least 1.4 kcal/mol relative to sites S₂ or S₄. 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 Ba²⁺ ions are bound simultaneously in sites S₂ and S₄. The x-ray structure was obtained from crystals that were grown from a solution with 5 mM BaCl₂, which is considerably larger than the micromolar Ba²⁺ concentrations used in the single-Ba²⁺ occupancy blockade experiments (Piasta et al., 2011). The implication is that single occupancy of site S₃ 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 Ba²⁺ in the binding sites is essential to definitively interpreting these experiments.

The simulations in this study were undertaken to help interpret the Ba²⁺ block experiments of the KcsA channel performed by Piasta et al. (2011). These calculations were designed to examine the translocation of the Ba²⁺ ion from site S₂ to site S₁ 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 Ba²⁺ ions in the Ba²⁺-soaked crystal structure of KcsA bind either to site S₂ or site S₄ of the selectivity filter. Second, the electrophysiological data suggested that the dominant situation during the prolonged blocks occurs when Ba²⁺ is bound to the outermost of these two sites. Under these conditions, the external lock-in site must either be site S₁ or site S₀, two possibilities that were investigated computationally.

Using umbrella sampling PMF calculations, we determined the affinity of K⁺ to sites S₀ and S₁ when site S₂ is occupied by either Ba²⁺ or K⁺. The computations indicate that binding of a cation (K⁺ or Na⁺) in site S₁ while a Ba²⁺ is bound in site S₂ is energetically prohibitive. Although site S₁ has a weak affinity for K⁺ when S₂ is occupied by K⁺, it is highly unfavorable for binding K⁺ when S₂ is occupied by Ba²⁺. This suggests that the proposal of Piasta et al. (2011) that the lock-in effect is due to a block of the translocation of Ba²⁺ from site S₂ to site S₁ caused by the binding of a K⁺ ion in the adjacent site S₁ is unlikely.

If S₁ is not the lock-in site, the only remaining external site that could serve in this capacity is S₀. This conclusion is somewhat unexpected and counterintuitive. Previous computational studies have indicated that S₀ is not selective for K⁺. In fact, this observation partly motivated Piasta et al. (2011) to rule out the possibility that S₀ 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 S₀ becomes considerably more selective when S₂ is occupied by Ba²⁺ than had been observed in previous studies when the filter was occupied by K⁺. The trend was confirmed by FEP calculations indicating that site S₀ sites became significantly more selective for K⁺ when Ba²⁺ 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 Ba²⁺ ion and the lock-in ion. This repulsion drives the lock-in ion bound to site S₀ 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 Ba²⁺ 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 Ba²⁺ permeation using Grote–Hynes theory and our calculated 1D PMF of Ba²⁺ translocation from S₂ to S₁. The calculated rate constant for this process is ∼52 s⁻¹, which is in reasonable agreement with the experimental rate of 204 s⁻¹. 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 Ba²⁺ permeates from S₂ to S₁ when there is an alkali ion in S₀. 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 Ba²⁺ permeation PMFs and estimated rate constant were consistent with this kind of lock-in scenario, the affinity of K⁺ to site S₀ 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 Ba²⁺ ion from site S₂ to site S₁, it is possible that the observed lock-in effects actually involve translocation steps with Ba²⁺ bound to sites that are more inward (e.g., S₃ or S₄). Ultimately, this will require a complete characterization of each translocation step involved when a Ba²⁺ 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.