TEA is a classical blocker of K+ channels. From mutagenesis studies, it has been shown that external blockade by TEA is strongly dependent upon the presence of aromatic residue at Shaker position 449 which is located near the extracellular entrance to the pore (Heginbotham, L., and R. MacKinnon. 1992. Neuron. 8:483–491). The data suggest that TEA interacts simultaneously with the aromatic residues of the four monomers. The determination of the 3-D structure of the KcsA channel using X-ray crystallography (Doyle, D.A., J.M. Cabral, R.A. Pfuetzner, A. Kuo, J.M. Gulbis, S.L. Cohen, B.T. Chait, and R. MacKinnon. 1998. Science. 280:69–77) has raised some issues that remain currently unresolved concerning the interpretation of these observations. In particular, the center of the Tyr82 side chains in KcsA (corresponding to position 449 in Shaker) forms a square of 11.8-Å side, a distance which is too large to allow simultaneous interactions of a TEA molecule with the four aromatic side chains. In this paper, the external blockade by TEA is explored by molecular dynamics simulations of an atomic model of KcsA in an explicit phospholipid bilayer with aqueous salt solution. It is observed, in qualitative accord with the experimental results, that TEA is stable when bound to the external side of the wild-type KcsA channel (with Tyr82), but is unstable when bound to a mutant channel in which the tyrosine residue has been substituted by a threonine. The free energy profile of TEA relative to the pore is calculated using umbrella sampling simulations to characterize quantitatively the extracellular blockade. It is found, in remarkable agreement with the experiment, that the TEA is more stably bound by 2.3 kcal/mol to the channel with four tyrosine residues. In the case of the wild-type KcsA channel, TEA (which has the shape of a flattened oblate spheroid) acts as an ideal plug blocking the pore. In contrast, it is considerably more off-centered and tilted in the case of the mutant channel. The enhanced stability conferred by the tyrosine residues does not arise from Π–cation interactions, but appears to be due to differences in the hydration structure of the TEA. Finally, it is shown that the experimentally observed voltage dependence of TEA block, which is traditionally interpreted in terms of the physical position of the TEA along the axis of the pore, must arise indirectly via coupling with the ions in the pore.

Introduction

An important fraction of the present information about the function of potassium channels has been deduced from studies with quaternary ammonium (QA) pore blockers, such as tetramethylammonium (TMA), TEA, tetrabutylammonium, and related alkyl derivatives (Armstrong 1968, Armstrong 1969, Armstrong 1971; French and Shoukimas 1981; Yellen 1998). In virtue of their ability to act as sensors of the pore of K+ channels, QAs are very important tools in electrophysiological studies of biological membranes.

All evidence support the idea that these agents act by blocking the permeation pathway, occluding the movement of K+ ions through the pore. Nearly all K+ channels are blocked by QAs on the intracellular side. The specificity of internal blockade is largely determined by the hydrophobicity of the blocker (Armstrong 1971; French and Shoukimas 1981). The observations are consistent with a blocker molecule that must move through part of the membrane voltage to reach its binding site inside the pore (Armstrong 1971; French and Shoukimas 1981). The mutations that most affect intracellular TEA blockade have been identified (MacKinnon and Yellen 1990; Yellen et al. 1991; Choi et al. 1993). In contrast, some K+ channels are blocked by TEA on the extracellular side, whereas others are not. The extracellular blockade is very selective for TEA over TMA. In addition, the affinity of extracellular TEA blockade is enhanced by the presence of an aromatic residue corresponding to position 449 in Shaker, which is presumably located near the extracellular mouth of the pore (MacKinnon and Yellen 1990; Heginbotham and MacKinnon 1992; Kavanaugh et al. 1992).

The recent determination of the 3-D structure of the KcsA channel from Streptomyces lividans using X-ray crystallography (Doyle et al. 1998) has raised some issues concerning the previous interpretation of electrophysiological data of extracellular blockade of K+ channels by TEA (Heginbotham and MacKinnon 1992). Those issues remain currently unresolved. In electrophysiological experiments, a linear relationship is observed between the free energy for TEA blockade and the number of channel subunits containing an aromatic residue at position 449, suggesting that these residues interact simultaneously with a TEA molecule (Heginbotham and MacKinnon 1992). Consistent with those results, it was proposed that TEA extracellular blockade was controlled by Π–cation interactions with a bracelet of pore-lining aromatic residues (see Fig. 7 of Heginbotham and MacKinnon 1992). However, this interpretation of electrophysiological data is somewhat at odds with the crystallographic structure of the KcsA channel (Doyle et al. 1998). As shown in Fig. 1, the center of the Tyr82 side chains in KcsA (which corresponds to position 449 in Shaker) of the four monomers form a square of 11.8-Å side. Because the diameter of TEA is roughly on the order of 6–8 Å, such a distance appears to be too large to allow simultaneous interactions of the TEA with the aromatic side chains. One assumes that KcsA, a bacterial channel, should be a relevant model. Based on a combination of structural and functional data with neurotoxin from scorpion, it has been concluded that KcsA is structurally similar to eukaryotic K+ channels such as Shaker (MacKinnon et al. 1998). Furthermore, like Shaker, KcsA is sensitive to external TEA blockade (Heginbotham et al. 1999; Meuser et al. 1999) and the affinity of external TEA blockade is decreased in Y82T KcsA mutant channels (Heginbotham et al. 1999). Nonetheless, the apparent inconsistency could be due to subtle structural differences between the pore of KcsA and Shaker. Alternatively, it is possible that considerations based upon a static view of the channel are insufficient. In particular, the influence of the aqueous environment and of the dynamics of the channel on TEA blockade may play an important role. Without further considerations, our current understanding of TEA extracellular blockade remains, therefore, incomplete.

Because of their widespread use in experimental studies of K+ channel, it is important to better understand the interaction and the association of QAs with K+ channels at the microscopic level. Theoretical studies based on molecular dynamics (MD) simulations of detailed atomic models provide a unique opportunity to understand the microscopic mechanism of QA blockades. The goal of this paper is to report results from MD simulations on extracellular TEA blockade of the KcsA K+ channel. In the last year or so, MD simulations of the KcsA channel with explicit solvent and membranes have been reported (Guidoni et al. 1999; Åqvist and Luzhkov 2000; Bernèche and Roux 2000; Luzhkov and Åqvist 2000; Shrivastava and Sansom 2000). Models of a wild-type KcsA (with Tyr at position 82) and of a mutant channel in which the residue has been substituted by a threonine were constructed, refined by energy minimization, and simulated with MD. Furthermore, the free energy profile of TEA relative to the entrance of the pore has been calculated using umbrella sampling simulations to quantitatively characterize the extracellular block. Our aim with these calculations is to extend our understanding of the microscopic mechanism by which TEA blocks K+ channels, and through this work, gain deeper insight into the structure and function of K+ channels.

Materials And Methods

Atomic Models and Simulation Method

All the MD simulations were performed using the program CHARMM (Brooks et al. 1983) with the all-atom PARAM22 potential function for proteins (MacKerell et al. 1998) and lipids (Schlenkrich et al. 1996) and the TIP3P model for water molecules (Jorgensen et al. 1983). The parameters of the positively charged TEA molecule were determined to be consistent with PARAM22. The atomic partial charges of TEA are as follows: N (−0.60), CA (−0.10), HA (0.25), CB (−0.27), and HB (0.09), where the first carbons (bound to the nitrogen) of TEA are labeled CA1–CA4, and the second carbons (bound to CA) are labeled CB1–CB4 (e.g., CB1 is bound to CA1). The hydrogens are labeled according to the carbon atom to which they are attached. The CA carbons are numbered in a cyclic order with the N-CA1 and the N-CA3 bonds pointing out of the plane and the N-CA2 and the N-CA4 bonds pointing into the plane. TEA can adopt two main conformations: (1) the fully symmetric (or quasi-planar) in which the CB1-CA1-N-CA3, CB2-CA2-N-CA4, CB3-CA3-N-CA1, and CB4-CA4-N-CA2 dihedral angles are all trans; and (2) the asymmetric (pyramidal) conformation in which the CB1-CA1-N-CA4, CB2-CA2-N-CA1, CB3-CA3-N-CA2, and CB4-CA4-N-CA3 dihedral angles are all trans. Ab initio quantum mechanical calculation of the intramolecular energy at the MP2/6–31G* level after geometry optimization at the HF/6–31G** level indicate that the symmetric conformer is more stable by 0.9 kcal/mol; the ab initio calculations were performed using the program Gaussian98 (Frisch et al. 1998). The electrostatic contribution to the solvation free energy was calculated using MD free energy simulations (MD/FES) of TEA in a sphere of 120 explicit water molecules with the spherical solvent boundary potential method (Beglov and Roux 1994). The charging procedure was done in the forward and backward directions with the windowing method: 11 windows with charge scaling were simulated with λi varying from 0.0 to 1.0 by an increment of 0.1 for a total of 220 ps. The result from the MD/FES windows was unbiased using the weighted histogram analysis method (Kumar et al. 1992). The calculated charging free energy of TEA is about −49 kcal/mol, in excellent accord with the experimental value of −49 kcal/mol (Aue et al. 1976) (the small positive free energy for creating the nonpolar cavity was not calculated). The calculated charging free energy difference between the symmetric and asymmetric conformer from MD/FES is ∼0.9 kcal/mol, favoring the symmetric conformer. Combining the results from the ab initio and solvation free energy calculations indicates that the asymmetric conformer is less stable by 1.8 kcal/mol. Because the total energy difference corresponds to populations of ∼5 and 95%, all the simulations with the KcsA channel were generated using the symmetric conformer. The symmetric conformer of TEA has the shape of a flattened oblate spheroid. When its “plane” is perpendicular to the channel axis, as shown in Fig. 1, TEA acts as an “ideal plug” blocking the extracellular pore entrance.

Atomic models of KcsA with TEA in the external pore entrance were constructed for a wild-type channel (KcsA-Y82) and a mutant channel in which Tyr82 is substituted by a Thr (KcsA-T82). The substituted residue is located just after the signature sequence common to K+ channels: TTVGYGDM-T(449) for Shaker, and TTVGYGDL-Y(82) for KcsA. The initial configuration of the system used in the present study was taken from an equilibrated configuration of the KcsA K+ channel embedded in a DPPC membrane bathed by a 150-mM KCl aqueous salt solution (Bernèche and Roux 2000). The general construction and equilibration protocol for this initial system has been described previously (Bernèche and Roux 2000). The Y82T mutant was constructed from the equilibrated wild-type system by mapping directly the coordinates of the corresponding atoms (backbone and Cβ) and building the remaining side chain atoms of Thr82 in a standard geometry using internal coordinates. The total number of atoms in the system is slightly above 40,000 (KcsA, 112 DPPC, 6,532 water molecules, 3 K+ in the pore, and 12 K+ and 23 Cl in the bulk solution). In both systems, the side chain of Glu71 is constructed in a protonated state to form a di-acid hydrogen bond with the carboxylate group of Asp80 near the extracellular surface (Roux et al. 2000). The entire system is electrically neutral. The channel axis is oriented along the z-axis; the center of the membrane is at z = 0.

TEA was first docked onto the extracellular entrance of the two channel models (KcsA-Y82 and KcsA-T82). The initial position of TEA was ∼21 Å along the z-axis (Fig. 1). The bulk water molecules overlapping with the TEA were removed. The crystallographic K+ ion in the outer site (located at z =17 Å) was replaced by a water molecule because it was <4 Å from the TEA. The crystallographic water and the K+ in the innermost binding site were kept. The initial configuration of the pore is W-K-W-W-TEA. In addition, there is one K+ in the water-filled cavity. After the initial construction, the atomic models were refined with energy minimization and equilibrated for 450 ps in the presence of progressively decreasing structural energy restraints to prevent any rapid and spurious displacements. All atomic restraints were gradually reduced to zero, such that the systems were completely free of restraints at the end of the equilibration (one center of mass planar harmonic restraint was applied along z to the bilayer during the production dynamics). Unrestrained trajectories of 1.75 and 1.2 ns were generated at 315 K for the KcsA-Y82 and KcsA-T82 systems, respectively. We refer to these two trajectories as Y82-MD1 and T82-MD1. To monitor the association of TEA with the external mouth of the channel, we calculated the instantaneous distance r between the center of mass of TEA and the center of mass of the selectivity filter (defined as the backbone of residues Thr75-Val76-Gly77-Tyr-78) as a function of time. The simulation of TEA in the Y82T mutated channel was stopped after the TEA escaped from the extracellular entrance of the channel and drifted away into the bulk solvent. The trajectories were generated at constant pressure with periodic boundary conditions (Feller et al. 1995). The electrostatic interactions were treated with no truncation using the particle mesh ewald method (Essmann et al. 1995).

The potential of mean force (PMF) was calculated using the “umbrella sampling method” (Valleau and Torrie 1977). This technique consists in simulating the system in the presence of an artificial biasing potential (the “umbrella potential”) to improve the sampling around some particular region of configurational space (Roux 1995); the calculations are then postprocessed to yield unbiased results. In the present case, umbrella sampling MD simulations were generated in the presence of a simple harmonic potential of the form:

\begin{equation*}w_{i} \left \left(r\right) \right =\frac{1}{2}K\; \left \left[r-r^{ \left \left(i\right) \right }\right] \right ^{2}{\mathrm{,}}\end{equation*}
1

where r is the distance between the center of mass of the TEA and the center of mass of the selectivity filter (backbone of residues Thr75-Val76-Gly77-Tyr-78). A total of 15 independent simulations, with a biasing harmonic potential successively centered on r(i) equal to 8.5, 9.0,…, and 15.5 Å, were generated. A force constant of K =10 kcal/mol·Å2 was used. The center-of-mass distance biasing potential was setup using the MMFP facility of CHARMM (Brooks et al. 1983). Each umbrella sampling trajectory was generated according to the same simulation protocol as described above, and the atomic coordinates of the system were saved for analysis. The initial system was reconstructed in the configuration K-W-K-W-TEA with two K+ ions in the selectivity filter on the basis of the results from the unrestrained simulations (see results and discussion). The initial models were refined with energy minimization and equilibrated for nearly 200 ps according to a similar protocol as described above. The windows were equilibrated and generated successively from 8.5 to 15.5 Å. The protocol for generating each window involved an equilibration period of 10 ps followed by production simulation and data collection. The PMF for the wild-type KcsA-Y82 channel, WY82(r), was calculated from MD simulations of 114 ps per window with an additional 54 ps for the central region (9.5 < r(i) < 12.5 Å) for a better convergence (total simulation of 168 ps). The PMF for the KcsA-T82 mutant channel, WT82(r), was calculated from MD simulations of 124 ps per window. The weighted histogram analysis method was used to unbias and recombine the results from the 15 individual biased umbrella sampling simulations (Kumar et al. 1992; Brooks and Nilsson 1993; Roux 1995).

Finally, to complement the information about the bound state, additional trajectories of 0.5 ns were generated starting the TEA in the free energy wells determined in the PMF for the KcsA-Y82 and KcsA-T82 channels without umbrella sampling biasing potential. We refer to these two trajectories as Y82-MD2 and T82-MD2. In addition, a reference 100-ps trajectory of TEA in a 20-Å sphere of 1,064 water molecules was generated with spherical solvent boundary potential to calculate the solvent radial distribution function.

All the calculations were done on the SGI Origin at the NCSA, on a cluster of Linux PC, and on a COMPAQ SC232 supercomputer at the CEA in Grenoble. 1 ns of simulation takes ∼200 h on the SGI Origin with 16 R10000 CPU running in parallel with the message passing interface.

Results And Discussion

The distance between the center of mass of TEA and the center of mass of the main residues forming the selectivity filter is shown as a function of time in Fig. 2 for the Y82-MD1 and T82-MD1 trajectories, respectively. It is observed that TEA remains bound to the extracellular mouth of the wild-type KcsA channel during the complete trajectory. In contrast, it escapes after 1.1 ns from the extracellular entrance of the Y82T mutant channel. Even before the escape event, TEA appears to be loosely bound to the mutant channel, whereas the association of TEA with the wild-type channel is very stable. The average distance r is 11.5 Å for the Y82-MD1 simulation, whereas it is 13.0 Å during the T82-MD1 trajectory. Therefore, TEA is nearly 3 Å deeper on average when it is bound to the channel with aromatic residues at position 82.

The overall channel conformation appears to be essentially unaffected by the amino acid substitution at position 82, though there are some small structural changes in the vicinity of the pore entrance. A superposition of the structures corresponding to averages over 500 ps of the two trajectories is shown in Fig. 3. The corresponding root-mean-squared (rms) difference between the backbone of the wild-type and mutant channels for residues 77–84 is on the order of 1.4 Å. In comparison, the deviation between the four different monomers of a given channel (Y82 or T82) varies between 0.3 and 0.9 Å for the backbone of residues 77–84. Therefore, the difference between the backbone of the wild-type and mutant channel is slightly larger than the normal deviations of the channel monomers from one another.

Even though the TEA molecule remains associated with the extracellular entrance of the KcsA-Y82 channel throughout the trajectory, it exhibits nonetheless significant dynamics, undergoing rapid rotational motions. In addition, a spontaneous concerted transition of the K+ ions in the inner site and in the water-filled cavity occurred within the first 100 ps of the production period during both the Y82-MD1 and T82-MD1 simulations. The initial configuration of the pore was W-K-W-W-TEA plus one K+ in the cavity (as described in materials and methods). During the transition, the K+ in the inner site translocated to the central site (which is occupied by a water molecule in the crystallographic structure), whereas the K+ which was in the cavity moved into the innermost site. After 100 ps of production dynamics, the selectivity filter is thus occupied by two K+ ions located in the innermost and central sites, respectively, whereas no ion is left in the water-filled cavity. The resulting configuration, K-W-K-W-TEA, is reminiscent of the configuration with three ions in the selectivity filter observed in previous MD simulations of the KcsA channel (Bernèche and Roux 2000).

The results from the trajectories described above are in qualitative accord with experimental observations of Heginbotham and MacKinnon 1992: extracellular TEA blockade is more stable in the presence of aromatic residues corresponding to Shaker position 449 near the channel entrance. Nonetheless, this type of analysis based on single trajectories has obvious limitations. To make progress, it is necessary to characterize extracellular TEA blockade quantitatively for the two channels. A powerful approach is to compute the free energy profile or the PMF, corresponding to the association of TEA with the extracellular entrance of the channel. The PMF, which was first introduced by Kirkwood 1935 in statistical mechanical theories of liquids, is a key concept in modern statistical mechanical theories of complex molecular systems (Brooks et al. 1988). In particular, binding constants can be can be conveniently expressed in terms of the PMF (Gilson et al. 1997; Roux and Simonson 1999). To compute the PMF, we performed umbrella sampling simulations (Valleau and Torrie 1977). These simulations are performed in the presence of the artificial biasing potential between the channel and the TEA molecule given in to enhance the configurational sampling. However, it must be stressed that the bias introduced by this potential is rigorously removed in postanalysis, thus, enabling us to characterize the unbiased free energy surface of the system.

The calculated PMF of TEA in the extracellular mouth of the KcsA-Y82 and KcsA-T82 channels, WY82(r) and WT82(r), are shown in Fig. 4 (r corresponds to the distance between the center of mass of TEA and the center of mass of the main residues forming the selectivity filter). For the KcsA-Y82 channel, there is a first free energy well at 9.8 Å followed by a second well ∼12 Å. For the KcsA-T82, there is a single free energy well located at 11.6 Å. In both cases, the TEA is bound more deeply into the channel mouth than was observed in the previous trajectories Y82-MD1 and T82-MD1 (see below). Although the two PMF calculated with umbrella sampling are determined up to an arbitrary offset constant (see Fig. 4 legend for more detail), a comparison of the relative free energy profile of TEA for the two channels is meaningful and the bound state is clearly more stable in the case of the KcsA-Y82. The relative binding free energy ΔGY82T = (GT82GY82) is given by:

\begin{equation*}{\mathrm{{\Delta}}}G_{{\mathrm{Y82T}}}=-k_{{\mathrm{B}}}T\;{\mathrm{ln}} \left \left[\frac{{\int _{0}^{r_{c}}}dr\;{\mathrm{exp}} \left \left({-W_{{\mathrm{T82}}} \left \left(r\right) \right }/{k_{{\mathrm{B}}}T}\right) \right }{{\int _{0}^{r_{c}}}dr\;{\mathrm{exp}} \left \left({-W_{{\mathrm{Y82}}} \left \left(r\right) \right }/{k_{{\mathrm{B}}}T}\right) \right }\right] \right {\mathrm{,}}\end{equation*}
2

where rc, chosen here as 14 Å, is the cutoff distance that defines the bound state. Evaluating the Boltzmann integrals in numerically using a trapezoidal rule yields a free energy difference of +2.3 kcal/mol. The variations in the PMF between 15 and 16 Å suggest that the uncertainty in the offset energy constant is on the order of 0.4 kcal/mol (Fig. 4). This provides an estimate of the uncertainty on the calculated ΔGY82T. The experimentally determined TEA/Shaker dissociation constant with four tyrosines at position 449 is 0.65 mM, whereas it is 22 mM after substitution by four threonines, corresponding to a loss of binding free energy of +2.1 kcal/mol (Heginbotham and MacKinnon 1992). Given the complexity of the molecular system being considered, the agreement with the experimental data is remarkable.

To better characterize the dynamical fluctuations of TEA in the free energy wells, we generated two additional trajectories of 0.5 ns for the two channels with no biasing potential (Y82-MD2 and T82-MD2). Typical conformations of TEA in the extracellular binding site are shown in Fig. 5 for the two channels. It is observed that in the case of the KcsA-Y82 channel, TEA acts as an ideal plug with its plane being nearly perpendicular to the channel axis and its end carbons (CB) pointing toward the tyrosine side chains. Analysis of the time course of the atoms shows that TEA maintains a given orientation for periods of a few hundred picoseconds before undergoing some rotational fluctuations. Because TEA is not perfectly spherical, the mode of interaction of TEA with the channel is affected by the rotation dynamics. The instantaneous position of the center of mass of TEA and of the Tyr82 side chains in the plane orthogonal to the channel axis are shown in Fig. 6. Although TEA remains well centered, on average, along the pore axis between the four tyrosine side chains, there are some lateral fluctuations on the order of 0.7 Å. The TEA molecule spends nearly an equal amount of time in contact with the residues of each of the four monomers. In contrast, it is considerably more tilted and off-centered during the T82-MD2 trajectory, spending >80% of its time near monomer A and the rest near monomer C (Fig. 5). Nonetheless, it is clear that the presence of TEA in the extracellular site of the mutant channel would prevent the passage of K+ ions despite the significant tilt and off-center position. As can be seen from Fig. 5, there is no room for the passage of ions between the TEA and the entrance to the pore.

The orientation adopted by TEA in the Y82-MD2 trajectory is consistent with experimental results with heterogenous tetrameric channels containing two tyrosines and two threonines, which suggest that the TEA molecule interact simultaneously with all the aromatic side chains (Heginbotham and MacKinnon 1992). Based on these observations, it was proposed that the increased affinity of TEA external blockade conferred by the aromatic residues Tyr82 arises from Π–cation interactions (Heginbotham and MacKinnon 1992). Such interactions, which govern the association of simple cations with benzene (Sunner et al. 1981; Gao et al. 1993), are also observed in proteins when a positively charged moiety is facing the plane of an aromatic side chain (Phe, Tyr, or Trp; Gallivan and Dougherty 1999). However, in the present case, the geometry is not optimal for Π–cation interactions: TEA does not face directly the plane of the aromatic side chain since the end carbons (CB) are closest to the CE2 carbons of the tyrosine ring (average distance of 4.5 Å with fluctuations on the order of 0.8 Å). In this configuration, the direct electrostatic interaction between the TEA molecule and each of the Tyr82 side chain fluctuate between positive and negative values, being roughly 0 kcal/mol on average withrms fluctuations of 2.3 kcal/mol. The corresponding van der Waals interactions are on the order of −0.8 kcal/mol. The unfavorable electrostatic interactions arises from the fact that the TEA molecule does not face the aromatic ring (Fig. 5), where the field from the electric quadrupole would be favorable. As a comparison, the most favorable interaction between TEA and a single isolated tyrosine side chain (without the backbone) is about −14 kcal/mol with the present potential function, a value which is similar to that of simple cations with benzene in the gas phase (Sunner et al. 1981). This interaction is achieved when TEA is facing directly the plane of the aromatic ring in an optimal configuration. The instantaneous electrostatic interactions observed between TEA and the Tyr82 side chains during the Y82-MD2 dynamical trajectories are much smaller than this value. Thus, this analysis shows that the increased affinity of TEA external blockade conferred by the aromatic residues Tyr82 relative to the Thr82 mutant does not arise directly from Π–cation interactions. What is then the microscopic origin of the enhanced affinity of TEA for the KcsA-Y82 channel?

Because the TEA is a large “hydrophobic” cation, hydration forces may play an important role in the stabilization of the extracellular blockade. To examine the hydration structure of TEA, we calculated the radial distribution of water oxygen around the nitrogen atom of TEA. The results are shown in Fig. 7. For comparison, the same distribution function calculated for TEA in bulk liquid water is also shown. It is clearly observed that the radial water density around TEA is highest for the bulk liquid environment, and lowest in the case of the KcsA-Y82 simulation. In the case of the KcsA-T82 simulation, the water density is somewhat intermediate between those of the bulk liquid and the wild-type channel. The number of water molecules in the first hydration shell (defined as 0 < r < 4.8 Å) is different: there are four and five water molecules around TEA bound to the KcsA-Y82 and KcsA-T82 channels, respectively. The number of water molecules in the second hydration shell (defined as 4.8 < r < 7.3 Å) differs sharply: there are 16 second-shell water molecules in the KcsA-Y82 simulation, whereas there are 22 in the KcsA-T82 simulation. For comparison, there are 8.5 water molecules in the first solvation shell of TEA in bulk water and 41 in the second shell. Thus, there is an excess of nearly seven water molecules around TEA bound to the mutant channel relative to the wild-type channel. The smaller hydration number is caused in part by the fact that TEA bound to the KcsA-Y82 channel lies like an ideal plug at the entrance of the pore in a pocket formed by the bulky tyrosine side chains. This leaves no extra space for water molecules. In contrast, TEA bound to the KcsA-T82 channel is considerably more tilted, and its hydration structure is not strongly altered by the relatively small threonine side chains. Variations in the hydration structure are known to be correlated with hydrophobic interaction between nonpolar molecules (Pangali et al. 1979a,Pangali et al. 1979b). This suggests that the hydration structure of TEA may be responsible for the difference in stability of extracellular blockade revealed by the PMF calculations.

According to the calculated PMF shown in Fig. 4, TEA is bound ∼1.5 Å deeper along the axis of the KcsA-Y82 channel than the KcsA-T82. Intriguingly, this result seems to contradict the increased voltage sensitivity of extracellular TEA blockade for Shaker that is observed experimentally when the aromatic residues at position 449 are substituted by threonines (Heginbotham and MacKinnon 1992); the weak voltage dependence increases from 4 to 19% upon the substitution. By interpreting the observed difference in terms of an “electric distance” (Hille 1992), it has been concluded that TEA binds less deeply into the transmembrane field when there are aromatic residues at position 449. However, further considerations indicate that this interpretation may be incorrect. The average position of the TEA along the axis of the pore is roughly ∼20 Å for the KcsA-Y82 and KcsA-T82 channels, respectively. Because of its size, it is impossible for the TEA to penetrate deeper than ∼19 Å into the narrow pore. But such a position is well outside the spatial region in which the transmembrane electric field is nonzero. It has been shown using continuum electrostatic calculations that the transmembrane potential varies in a nonlinear fashion between −18 and +18 Å, with the largest fraction of the potential dropping across the selectivity filter (Roux et al. 2000). The transmembrane field is vanishingly small beyond 18 Å. Therefore, there cannot be any direct interaction between the transmembrane voltage and the charged TEA when it is bound in the extracellular sites. On this basis, we conclude that the experimentally observed voltage dependence of TEA block must arise via a coupling with the ions in the pore. This result is intriguing and suggests many questions about the origin of the observed voltage dependence and the nature of the interaction of TEA with the ions in various occupancy state of the channel. For example, further analysis reveals that there is only one water molecule between the TEA and the nearest K+ ion in the selectivity filter (K-W-K-W-TEA) for the configurations corresponding to the first free energy well, whereas there can be one or two water molecules for the second well (K-W-K-W-W-TEA). In contrast, there is always only a single water molecule between TEA and the nearest K+ ion (K-W-K-W-TEA) for the KcsA-T82 channel. How these water molecules shield the repulsion from the K+ in the filter is not known. Undoubtedly, the complex coupling between the K+ located inside the selectivity filter will require further investigations to draw a complete picture of external blockade by TEA.

One may note that the TEA is clearly bound more stably during both of the MD2 trajectories than it was in the MD1 trajectories, which were described above (Fig. 2 and Fig. 4). In the case of the KcsA-Y82 channel, the average position of TEA is nearly 1–1.5 Å deeper in the MD2 simulation than in the MD1 simulations. Obviously, the results are influenced by the fact that those trajectories were started with the TEA molecule in a very stable position, i.e., effectively at the bottom of a free energy well as determined by the PMF calculations (Fig. 4). In contrast, the initial simulations MD1 were started in a plausible conformation, without any prior knowledge of the PMF. The initial starting position were, thus, imperfect. However, it should be stressed that, even though these initial simulations were limited, they nevertheless provided essential information in our analysis of extracellular TEA blockade. This highlights the importance of following up preliminary explorations based on unbiased simulations with quantitative analysis based on free energy PMF calculations to reach robust conclusions. Therefore, the present study provides a good illustration of the type of computational strategy that is required to make progress in simulation studies of biological macromolecular system.

After the present article was submitted, a study of TEA block was published by Luzhkov and Åqvist 2001. The binding of TEA in the symmetric (quasi-planar) and asymmetric form (pyramidal) was examined using automated docking and simulations of a reduced model. In accord with the present results, the origin of the increased stabilization by aromatic side chain was attributed to hydrophobic interactions. However, in sharp disagreement with the present results, it was concluded that the association of the symmetric conformer of TEA with the channel is unstable, and that external blockade is dominated by the binding of the asymmetric form of TEA. This result is puzzling since high level ab initio calculations and MD/FES indicate that the asymmetric form is poorly populated in solution (materials and methods). Luzhkov and Åqvist 2001 based their conclusions on free energies estimated by the linear interaction energy approximation. The present conclusions are based on the free energy profiles shown in Fig. 4. They were calculated using umbrella sampling simulations, which is an exact sampling method (Valleau and Torrie 1977). It seems likely that this is the origin of the discrepancy.

Conclusion

We have performed MD simulations of external blockade by TEA for two K+ channels: the wild-type KcsA, and a single point mutant channel in which Tyr82 has been substituted by a threonine. The simulations show that the stability of external blockade by TEA observed in the KcsA-Y82 and KcsA-T82 simulations is in qualitative agreement with experimental results: increased stability of external blockade is correlated with the presence of an aromatic residue at position 82 (corresponding to 449 in Shaker). To go beyond qualitative observations and characterize the extracellular blockade quantitatively, we calculated the free energy profile of TEA along the channel axis using umbrella sampling simulations. According to the PMF calculations, the free energy loss upon substitution of the four tyrosines by threonines is ∼2.3 ± 0.4 kcal/mol, in remarkable agreement with the experimental results of Heginbotham and MacKinnon 1992. Although this was not immediately obvious from considerations based only on a static view, these calculations show that the crystallographic structure of the KcsA K+ channel is consistent with the experimentally observed influence of aromatic side chains on extracellular TEA blockade. This further demonstrates the structural similarity of the extracellular entrance to the pore of KcsA with that of Shaker, thus, reinforcing the KcsA channel as a relevant model to help understand a large class of biologically important channels.

Most dramatically, the dynamical simulations provide a novel perspective of TEA blockades of K+ channels that challenges current views based on static configurations. In particular, although TEA remains well centered along the pore axis on average, there are significant lateral fluctuations, and it is not inserted tightly between the aromatic Tyr82 side chains of the four monomers. As a result, it does not simultaneously make Π–cation interactions with those residues as previously proposed (Heginbotham and MacKinnon 1992). When there are four tyrosines in position 82, the TEA maintains an average orientation that resembles that of a plug blocking the pore (Fig. 5). In contrast, the TEA is fairly tilted and more outwardly bound when the four tyrosines are substituted by threonines. Although this appears to contradict the interpretation of the experimentally observed voltage dependence of TEA block in terms of the physical position of the TEA along the axis of the pore, further considerations show that the voltage dependence must arise indirectly via a coupling with the ions in the pore. Clearly, further investigations will be required to clarify the interactions between the TEA and the K+ located inside the selectivity filter and the coupling to the transmembrane voltage.

The present results show that computations based on detailed atomic models can contribute to our understanding of the microscopic nature of extracellular blockade of K+ channels by TEA. Future work will be used to investigate extracellular blockade by other common QAs such as TMA and tetrabutylammonium and the influence of other mutations. Furthermore, the relative stability of TEA bound to various heterotetramer channels, with two tyrosines and two threonines residues in different subunits, will be examined and compared with experiments (Heginbotham and MacKinnon 1992). Finally, the free energy profile will be extended to large distances to calculate the absolute dissociation constant of TEA in the extracellular binding site.

Acknowledgments

We are grateful to R. MacKinnon and L. Heginbotham for helpful discussions, and to A.D. MacKerell for his help with the ab initio calculations.

This work was supported by the National Institutes of Health grant R01-GM62342-01(to B. Roux).

References

Åqvist
J.
,
Luzhkov
V.
Ion permeation mechanism of the potassium channel
Nature.
404
2000
881
884
[PubMed]
Armstrong
C.M.
Induced inactivation of the potassium permeability of squid axon membranes
Nature.
219
1968
1262
1263
[PubMed]
Armstrong
C.M.
Inactivation of the potassium conductance and related phenomena caused by quaternary ammonium ion injection in squid axons
J. Gen. Physiol.
54
1969
553
575
[PubMed]
Armstrong
C.M.
Interaction of tetraethylammonium ion derivatives with the potassium channels of giant axons
J. Gen. Physiol.
58
1971
413
437
[PubMed]
Aue
D.H.
,
Webb
H.M.
,
Bowers
M.T.
A thermodynamic analysis of solvation effects on the basicity of alkylamine. An electrostatic analysis of subtituents effects
J. Am. Chem. Soc.
98
1976
318
329
Beglov
D.
,
Roux
B.
Finite representation of an infinite bulk systemsolvent boundary potential for computer simulations
J. Chem. Physiol.
100
1994
9050
9063
Bernèche
S.
,
Roux
B.
Molecular dynamics of the KcsA K(+) channel in a bilayer membrane
Biophys. J.
78
2000
2900
2917
[PubMed]
Brooks
B.R.
,
Bruccoleri
R.E.
,
Olafson
B.D.
,
States
D.J.
,
Swaminathan
S.
,
Karplus
M.
CHARMMa program for macromolecular energy minimization and dynamics calculations
J. Comp. Chem.
4
1983
187
217
Brooks
C.L.
III
,
Nilsson
L.
Promotion of helix formation in peptides dissolved in alcohol and water-alcohol mixtures
J. Am. Chem. Soc.
115
1993
11034
11035
Brooks, C.L., III, M. Karplus, and B.M. Pettitt. 1988. Proteins. A theoretical perspective of dynamics, structure and thermodynamics. In Advances in Chemical Physics. I. Prigogine and S.A. Rice, editors. Vol. LXXI. John Wiley & Sons, New York.
Choi
K.L.
,
Mossman
C.
,
Aube
J.
,
Yellen
G.
The internal quaternary ammonium receptor site of Shaker potassium channels
Neuron.
10
1993
533
541
[PubMed]
Doyle
D.A.
,
Cabral
J.M.
,
Pfuetzner
R.A.
,
Kuo
A.
,
Gulbis
J.M.
,
Cohen
S.L.
,
Chait
B.T.
,
MacKinnon
R.
The structure of the potassium channelmolecular basis of K+ conduction and selectivity
Science.
280
1998
69
77
[PubMed]
Essmann
U.
,
Perera
L.
,
Berkowitz
M.L.
,
Darden
T.
,
Lee
H.
,
Pedersen
L.G.
A smooth particle mesh Ewald method
J. Chem. Phys.
103
1995
8577
8593
Feller
S.E.
,
Zhang
Y.H.
,
Pastor
R.W.
,
Brooks
B.R.
Constant pressure molecular dynamics simulationthe Langevin piston method
J. Chem. Physiol.
103
1995
4613
4621
French
R.J.
,
Shoukimas
J.J.
Blockage of squid axon potassium conductance by internal tetra-N-alkylammonium ions of various sizes
Biophys. J.
34
1981
271
291
[PubMed]
Frisch
M.J.
,
Trucks
G.W.
,
Schlegel
H.B.
,
Scuseria
G.E.
,
Robb
M.A.
,
Cheeseman
J.R.
,
Zakrzewski
V.G.
,
Montgomery
J.A.
,
Stratmann
R.E.
Gaussian 98
1998
Gaussian Inc
Pittsburgh, PA
Gallivan
J.P.
,
Dougherty
D.A.
Cation-pi interactions in structural biology
Proc. Natl. Acad. Sci. USA.
96
1999
9459
9464
[PubMed]
Gao
J.
,
Chou
L.W.
,
Auerbach
A.
The nature of cation-pi bindinginteractions between tetramethylammonium ion and benzene in aqueous solution
Biophys. J.
65
1993
43
47
[PubMed]
Gilson
M.K.
,
Given
J.A.
,
Bush
B.L.
,
McCammon
J.A.
The statistical-thermodynamic basis for computation of binding affinitiesA critical review
Biophys. J.
72
1997
1047
1069
[PubMed]
Guidoni
L.
,
Torre
V.
,
Carloni
P.
Potassium and sodium binding to the outer mouth of the K+ channel
Biochemistry.
38
1999
8599
8604
[PubMed]
Heginbotham
L.
,
MacKinnon
R.
The aromatic binding site for tetraethylammonium ion on potassium channels
Neuron.
8
1992
483
491
[PubMed]
Heginbotham
L.
,
LeMasurier
M.
,
Kolmakova-Partensky
L.
,
Miller
C.
Single Streptomyces lividans K(+) channels. Functional asymmetries and sidedness of proton activation
J. Gen. Physiol.
114
1999
551
560
[PubMed]
Hille
B.
Ionic Channels of Excitable Membranes
2nd ed
1992
Sinauer Associates, Inc
Sunderland MA
pp. 607
Jorgensen
W.L.
,
Chandrasekhar
J.
,
Madura
J.D.
,
Impey
R.W.
,
Klein
M.L.
Comparison of simple potential functions for simulating liquid water
J. Chem. Physiol.
79
1983
926
935
Kavanaugh
M.P.
,
Hurst
R.S.
,
Yakel
J.
,
Varnum
M.D.
,
Adelman
J.P.
,
North
R.A.
Multiple subunits of a voltage-dependent potassium channel contribute to the binding site for tetraethylammonium
Neuron.
8
1992
493
497
[PubMed]
Kirkwood
J.G.
Statistical mechanics of fluid mixtures
J. Chem. Physiol.
3
1935
300
313
Kumar
S.
,
Bouzida
D.
,
Swendsen
R.H.
,
Kollman
P.A.
,
Rosenberg
J.M.
The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method
J. Comp. Chem.
13
1992
1011
1021
Luzhkov
V.B.
,
Åqvist
J.
A computational study of ion binding and protonation states in the KcsA potassium channel
Biochim. Biophys. Acta.
1481
2000
360
370
[PubMed]
Luzhkov
V.B.
,
Åqvist
J.
Mechanisms of tetraethylammonium ion block in the KcsA potassium channel
FEBS Lett.
495
2001
191
196
[PubMed]
MacKerell
A.D.
Jr.
,
Bashford
D.
,
Bellot
M.
,
Dunbrack
R.L.
,
Evanseck
J.D.
,
Field
M.J.
,
Fischer
S.
,
Gao
J.
,
Guo
H.
,
Joseph-McCarthy
D.
All-atom empirical potential for molecular modeling and dynamics studies of proteins
J. Physiol. Chem. B.
102
1998
3586
3616
MacKinnon
R.
,
Yellen
G.
Mutations affecting TEA blockade and ion permeation in voltage-activated K+ channels
Science.
250
1990
276
279
[PubMed]
MacKinnon
R.
,
Cohen
S.L.
,
Kuo
A.
,
Lee
A.
,
Chait
B.T.
Structural conservation in prokaryotic and eukaryotic potassium channels
Science.
280
1998
106
109
[PubMed]
Meuser
D.
,
Splitt
H.
,
Wagner
R.
,
Schrempf
H.
Exploring the open pore of the potassium channel from Streptomyces lividans
FEBS Lett.
462
1999
447
452
[PubMed]
Pangali
C.
,
Rao
M.
,
Berne
B.J.
A Monte Carlo simulation of the hydrophobic interaction
J. Chem. Physiol.
71
1979
2975
2981
a
Pangali
C.
,
Rao
M.
,
Berne
B.J.
Hydrophobic hydration around a pair of apolar species in water
J. Chem. Physiol.
71
1979
2982
2990
b
Roux
B.
The calculation of the potential of mean force using computer simulations
Comp. Physiol. Commun.
91
1995
275
282
Roux
B.
,
Simonson
T.
Implicit solvent models
Biophys. Chem.
78
1999
1
20
[PubMed]
Roux
B.
,
Bernèche
S.
,
Im
W.
Ion channels, permeation and electrostaticsinsight into the function of KcsA
Biochemistry.
39
2000
13295
13306
[PubMed]
Schlenkrich, M.J., J. Brickmann, A.D. MacKerell, Jr., and M. Karplus. 1996. An empirical potential energy function for phospholipids: criteria for parameters optimization and applications. In Biological Membranes. A Molecular Perspective from Computation and Experiment. K.M. Merz and B. Roux, editors. Birkhauser, Boston. 31–81.
Shrivastava
I.H.
,
Sansom
M.S.P.
Simulation of ion permeation through a K channelmolecular dynamics of KcsA in a phospholipid bilayer
Biophys. J
78
2000
557
570
[PubMed]
Sunner
J.
,
Nishizawa
K.
,
Kebarle
P.
Ion-solvent molecule interactions in the gas phase. The potassium ion and benzene
J. Physiol. Chem.
85
1981
1814
1820
Valleau
J.P.
,
Torrie
G.M.
,
A guide for monte carlo for statistical mechanics
Berne
B.J.
Statistical Mechanics, Part A
1977
169
194
Plenum Press
NY
Yellen
G.
The moving parts of voltage-gated ion channels
Q. Rev. Biophys.
31
1998
239
295
[PubMed]
Yellen
G.
,
Jurman
M.E.
,
Abramson
T.
,
MacKinnon
R.
Mutations affecting internal TEA blockade identify the probable pore-forming region of a K+ channel
Science.
251
1991
939
942
[PubMed]

Abbreviations used in this paper: MD, molecular dynamics; PMF, potential of mean force; QA, quaternary ammonium; rms, root-mean-squared; TMA, tetramethylammonium.