Article Information

• PDF (424 kb)
• Full Text with Large Figures
• Cited by in Scopus (1)

PubMed

• Articles by Matthew Hoyles
• Articles by Shin-Ho Chung

Related Articles

• Glycerol Conductance and Physical Asymmetry of the Escherich...

  Glycerol Conductance and Physical Asymmetry of the Escherichia coli Glycerol Facilitator GlpF Biophysical Journal, Volume 85, Issue 5, 1 November 2003, Pages 2977-2987 Deyu Lu, Paul Grayson and Klaus Schulten Abstract The aquaglyceroporin GlpF is a transmembrane channel of that facilitates the uptake of glycerol by the cell. Its high glycerol uptake rate is crucial for the cell to survive in very low glycerol concentrations. Although GlpF allows both influx and outflux of glycerol, its structure, similar to the structure of maltoporin, exhibits a significant degree of asymmetry. The potential of mean force characterizing glycerol in the channel shows a corresponding asymmetry with an attractive vestibule only at the periplasmic side. In this study, we analyze the potential of mean force, showing that a simplified six-step model captures the kinetics and yields a glycerol conduction rate that agrees well with observation. The vestibule improves the conduction rate by 40% and 75% at 10-M and 10-mM periplasmic glycerol concentrations, respectively. In addition, neither the conduction rate nor the conduction probability for a single glycerol (efficiency) depends on the orientation of GlpF. GlpF appears to conduct equally well in both directions under physiological conditions. Abstract | Full Text | PDF (407 kb)

• Computational Studies of Membrane Channels

  Computational Studies of Membrane Channels Structure, Volume 12, Issue 8, 1 August 2004, Pages 1343-1351 Benoit Roux and Klaus Schulten Summary The determination of the structure of several members of the K channel and aquaporin family represents a unique opportunity to explain the mechanism of these biomolecular systems. With their ability to go beyond static structures, molecular dynamics simulations offer a unique route for relating functional properties to membrane channel structure. The recent progress in this area is reviewed. Summary | Full Text | PDF (270 kb)

• A Gate in the Selectivity Filter of Potassium Channels

  A Gate in the Selectivity Filter of Potassium Channels Structure, Volume 13, Issue 4, 1 April 2005, Pages 591-600 Simon Bernèche and Benoît Roux Summary The selectivity filter of potassium channels is the structural element directly responsible for the selective and rapid conduction of K, whereas other parts of the protein are thought to function as a molecular gate that either permits or blocks the passage of ions. However, whether the selectivity filter itself also possesses the ability to play the role of a gate is an unresolved question. Using free energy molecular dynamics simulations, it is shown that the reorientation of two peptide linkages in the selectivity filter of the KcsA K channel can lead to a stable nonconducting conformational state. Two microscopic factors influence the transition toward such a conformational state: the occupancy of one specific cation binding site in the selectivity filter (S), and the strength of intersubunit interactions involving the GYG signature sequence. These results suggest that such conformational transitions occurring in the selectivity filter might be related to different K channel gating events, including C-type (slow) inactivation. Summary | Full Text | PDF (348 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 94, Issue 2, 366-378, 15 January 2008

doi:10.1529/biophysj.107.115139

Biophysical Theory and Modeling

Brownian Dynamics Theory for Predicting Internal and External Blockages of Tetraethylammonium in the KcsA Potassium Channel

Matthew Hoyles^*, Vikram Krishnamurthy^†, May Siksik^† and Shin-Ho Chung^*, ,  

^* Research School of Biological Sciences, The Australian National University, Canberra, Australia
^† Department of Electrical Engineering, The University of British Columbia, Vancouver, British Columbia, Canada

Address reprint requests to Shin-Ho Chung.

Tetraethylammonium (TEA) has long known to be a potent blocker of potassium currents flowing across the cell membrane. TEA has been widely used to separate different components of the membrane currents and to probe the dynamics of ion permeation across potassium channels. Many potassium channels are blocked by TEA molecules applied to the internal 1 and external sides of the membrane 2,3. Different K⁺ channels display different sensitivity to external TEA. For example, the inhibitory constant k_i for the Shaker K⁺ channel is 27mM, but its TEA affinity increases drastically when the threonine residue at position 449 is substituted with the tyrosine or phenylalanine residue 3,4. In KcsA and Kv2.1, where four tyrosine residues (at position 82 or 380), located just external to the selectivity filter, form a binding pocket for TEA, the current is halved when its concentration reaches between 1.7 and 4.5mM 5,6,7,8. The external TEA affinity in KcsA is drastically reduced when tyrosine at position 82 is substituted with nonaromatic residues, such as cysteine or threonine 5. The precise mechanism that influences TEA potency in different channels and the structural basis of the forces that stabilize a TEA molecule in the outer vestibule are not fully understood.

The one-dimensional free energy profile computed using molecular dynamics shows that a TEA molecule sees a shallow attractive potential well as it enters the binding pocket 9,10. This well extends ∼6Å in the direction normal to the membrane, and its deepest position, located outside of the selectivity filter, is ∼4.45kcal/mol or 7.5kT (1 kcal/mol equals 1.69kT at 298°C) 9. From steric considerations, it is postulated that the stabilizing effects between the aromatic rings of the tyrosine residues and a TEA molecule arise from hydration rather than electrostatic forces 9,10. On the other hand, using the nonsense suppression method in the mutant Shaker, Ahern et al. 11 provide experimental and computational evidence supporting the idea that external TEA blockade in Shaker arises, in large part, from a cation-π electron interaction between the aromatic ring of the tyrosine residues at position 449 and a quaternary ammonium ion.

The picture emerging from these studies is that in the KcsA and Shaker potassium channels, a TEA molecule physically occludes the pore when it is sequestered by four tyrosine residues located just external to the selectivity filter. In the Kv2.1 channel, which is a slowly inactivating delayed rectifier responsible for repolarization of the action potential, evidence suggests that current block by external TEA may hinge on several other factors. For example, TEA is ineffective in blocking Na⁺ currents flowing across this channel, in the absence of K⁺6. The substitution of the tyrosine residues at position 380 (corresponding to position 449 in Shaker) with cysteines has relatively small effect on TEA block; this mutation shifts the inhibitory constant from 3 to 9mM 12. Thus, the stabilization of TEA in the external vestibule is influenced by, among other factors, the conformation of the vestibule, the electrostatic milieu near the binding pocket, and the occupancy state of the selectivity filter. It is possible that, in the course of a prolonged activation of the channel, TEA efficacy undergoes dynamic changes, as does the conformation of the external vestibule 13,14 and of the selectivity filter 15.

The aim of this article is to show that Brownian dynamics (BD) simulations can reveal further insights into the action of internally and externally applied TEA on the permeation dynamics. The BD simulations are run for a period long enough to study the effects of TEA in the simulation assembly on the conductance properties of the KcsA K⁺ channel. The ability to run simulations over a period long enough to determine statistically consistent estimates of the current flow across ion channels is a distinct advantage of BD over molecular dynamics. This is achieved by making two approximations. First, BD assumes a fixed protein structure, so that the explicit dynamics of protein molecules are not simulated. Second, in BD, the dynamics of water molecules are not simulated. Instead, the net effect of incessant collisions between water molecules and ions, or water molecules and TEA molecules, are represented as frictional and random forces. These random forces are modeled as Brownian motion based on a functional central limit theorem approximation. This means that the hydration structure of molecules in the simulation assembly, or the attractive interaction between molecules via hydration forces, are not taken into account.

The forces arising from the interaction between a charged TEA molecule and the quadrupolar moment of the aromatic ring are computed from the three-dimensional free energy profile a TEA molecule encounters as it moves from the outer vestibule toward the selectivity filter of the KcsA channel. As described in Appendix 1 , the one-dimensional free energy profile, taken from molecular dynamics calculations 9, is converted into a one-dimensional potential energy profile to be used in BD simulations. The shape of the one-dimensional potential energy profile for TEA remains virtually identical to that of the one-dimensional free energy profile but its depth increases slightly once these conversions are made.

With this modification of the free energy profile, we carry out BD simulations of the KcsA potassium channel with KCl ions and TEA molecules both present. In this novel BD algorithm, the motion of TEA is influenced by the electrostatic, random, and frictional forces, as well as the attractive force arising from its interaction with four tyrosine residues just outside of the selectivity filter. This free energy profile, invisible to other charged particles in the simulation assembly, is in turn modulated by the locations of resident K⁺ ions in the selectivity filter. With these simulation techniques, we uncover several additional insights into the mechanisms of internal and external TEA blockade. First, we reproduce the macroscopically observed binding affinity of TEA through the fundamental processes operating in electrolyte solutions. By constructing the current-TEA-concentration curves using different energy profiles with different well-depths, we are able to determine the binding energy between TEA and the channel to account for the experimentally observed inhibitory constant k_i. Second, we show how the inhibitory constant is pronouncedly influenced by the applied membrane potential and ionic concentrations of K⁺ ions in the reservoirs. Such dynamic interactions arising from the presence of many particles cannot be revealed from a static profile of potential of mean force. Third, with the BD simulation technique, we have been able to investigate the kinetics of TEA binding and unbinding. We show that they follow a first-order Markovian process, and construct the transition probability matrices at various applied potentials and ligand concentrations. The opening and closing of a single channel is governed by a Markovian process 16,17. Within each brief open episode in the presence of TEA, the channel is blocked and unblocked at a fast timescale, according to another Markovian process. Finally, we show how the channel is blocked by TEA applied to the internal side of the membrane.

The channel structure we use for this study is the same as that used in previous BD simulations of the KcsA potassium channel 18. The model channel is based on the crystal structure of Doyle et al. 19 with the pore expanded at the internal end of the channel using a cylindrical repulsive potential in molecular dynamics to have a radius of 5Å. As described previously, the selectivity filter is also expanded slightly such that the minimum radius is 1.4Å. The full charges of −e and +e are assigned to D₈₀ and R₆₄ located at the extracellular aspect of the channel, whereas a partial charge of −0.5e and +0.5e is assigned to the ionized pair, E₁₁₈ and R₁₁₇, guarding the intracellular gate. One unpaired arginine residue (Arg₂₇) and the two pairs of ionizable residues, namely, Glu₅₁/Arg₅₂ and Glu₇₁/Arg₈₉, are kept neutral, for the reasons detailed in Chung et al. 23. For this study, we selected this model, among several other permeable forms of KcsA (see, for example, 18), as it more accurately replicates the experimental data, such as the observed current-voltage-concentration profiles.

The algorithm for performing BD simulations is conceptually simple. The position and velocity of each individual ion in and around the ion channel evolves according to a continuous time stochastic dynamical system. The position and velocity of the i^th ion at time t with mass m and charge q located at a given position are determined by the force acting on it and are computed by integrating the equation of motion, known as the Langevin equation, a stochastic version of Newton’s equation of motion. Thus,

(1)

(2)

(3)

Ions are initially assigned random positions with a specified mean concentration and Boltzmann-distributed velocities within cylindrical reservoirs of 30Å radius attached to either end of the channel that mimic the intracellular and extracellular space. We calculate the total force acting on each ion in the assembly and then determine new positions for the ions a short time later using an algorithm described elsewhere 20,21,22. A multiple time-step algorithm is used, where a time-step of Δt=100 fs is employed in the reservoirs and 2 fs in the channel where the forces change more rapidly.

The electrostatic forces F_λ acting upon the ions are calculated by solving Poisson’s equation using a boundary element method 24. To do this we assign a dielectric constant of ϵ_p=2 to the protein and membrane and ϵ_w=60 to the water inside the channel as in previous BD studies of the KcsA channel 18. Since calculating the electric forces at every step in the simulation is very time-consuming, we store precalculated electric fields and potentials due to one- and two-ion configurations in a system of lookup tables and interpolate values for the given ion positions from these during the simulation 21.

For this purpose, the total electric potential ϕ_i encountered by a charged particle i (including TEA molecules) is decomposed into four components:

(4)

The short-range forces are used to keep the ions and TEA molecules in the simulation assembly and also to mimic other interactions between charged particles that are not included in the Coulomb interaction. To prevent charged particles from leaving the system, a hard-wall potential is activated when they are within one ionic radius of the reservoir boundaries. Thus, a particle colliding with the boundary of the reservoir is elastically scattered. For the ion- or TEA-wall interaction U^IW, we use the 1/r⁹ repulsive potential

(5)

Ideally, each TEA molecule should be represented as a rigid oblate shape, a spheroid generated by the revolution of an ellipse, with the long and short axes equal to 4.6 and 2.3Å. For this study, we use a simple, spherical model of TEA with the radius of either 2.3Å (for externally applied TEA) or 3.4Å (for internally applied TEA). The molecular weight and diffusion coefficients used for the spherical TEA model are, respectively, 130.25 and 0.87×10⁻⁹ m² s⁻¹29.

Crouzy et al. 9 calculated a free energy profile encountered by a TEA molecule as it moves from the external reservoir toward the selectivity filter (Fig. 1, top). To convert this one-dimensional free energy profile, calculated using CHARMM, to a three-dimensional free energy profile that can be incorporated into the BD algorithm, we first fit an analytical function to the profile obtained by Crouzy et al. 9 (Fig. 1, bottom, open circles). The digitized points are fitted with the sum of three simple algebraic functions:

(6)

(7)

(8)

Display large version of this figure

Figure 1

Model potassium channel. (Top) Two of the four subunits of the full experimentally determined KcsA protein are represented as ribbons. Water molecules are shown in gold. The circled region of the protein, containing the selectivity filter is enlarged on the right-hand side and shown in a stick-and-ball model, together with three water molecules and 2K⁺ ions (gold). In simulations, each TEA molecule in the external reservoir is represented as a sphere with radius of 2.3Å, and that in the internal reservoir as a sphere with radius of 3.2Å. (Bottom) The one-dimensional free energy profile obtained by Crouzy et al. 9 replotted (open circles). The digitized data points are fitted with three analytical functions (solid line). The origin of the profile is at the center of the mass of the selectivity filter.

We then transform the analytical one-dimensional free energy profile G_1d into a three-dimensional free energy profile G_3d, taking into account the effective cross-sectional area at each z position. This conversion is done by first obtaining a one-dimensional potential profile U^PMF and then adding to it the short-range potential which represents repulsion from the channel wall. The relationship between the free energy in one-dimension G_1d and that in three-dimension G_3d is given in Appendix 1 . The depth of the well of the one-dimensional potential profile U^PMF is ∼2.5 kT deeper than the free energy profile in one-dimension G_1d.

The position of the center of mass of the selectivity filter, which is taken as the origin of the z axis in Fig. 1, is at z=15.1Å in our coordinate system (where z=0 is taken to be at the midway between the two ends of the channel). The three-dimensional free energy profile G_3d is spliced onto the usual BD force field using a reference potential U^REF(z). This potential is used along the z axis in the region (z=24.1 to 31.1Å) with two K⁺ ions in the selectivity filter of the channel to reproduce as closely as possible the conditions of the MD simulation. Except for short-range potentials, all electrostatic effects, such as direct and indirect image effects of the channel dielectric boundary, are included. The reference potential is then subtracted to avoid double counting of potentials already included from the molecular dynamics free energy profile. The positions of the two K⁺ ions in the selectivity filter are allowed to vary to minimize the total energy. These two positions range from z=9.55 to 9.71Å, and from z=16.47 to 16.66Å. This reference potential is then subtracted from the BD potential where the potential profile U^PMF is added. When the configuration of ions is different, these differences are taken into account by the resulting change in the BD potential.

We employ, when needed to maintain low concentrations of TEA in the external reservoir, the grand canonical Monte Carlo method 31. The effective volume of each reservoir in our simulation system is 8.7×10⁻²⁶m³. Since placing one molecule in the reservoir will bring its concentration to ∼20mM, we make use of this method to maintain a low TEA concentration, rather than increasing the size of the reservoirs. The technical details of implementing the grand canonical Monte Carlo method are given in Corry et al. 32.

Our BD simulations reveal that internally applied TEA acts by physically occluding the pore. BD simulations are carried out with either one or more TEA molecules in the internal vestibule, and 16K⁺ and 16 Cl⁻ ions in each reservoir. In these simulations, the free energy profile calculated from molecular dynamics is not included. Yellen et al. 2 identified an amino-acid residue in the Drosophila Shaker channel that influences the affinity for TEA applied intracellularly. The current in the wild-type Shaker (ShIR) is halved when 1mM TEA is applied on the internal face of the channel. Substitution of threonine residue at position 441 with serine drastically alters sensitivity to internal TEA. The sensitivity of the T441S mutant channel to TEA is reduced by an order of magnitude compared to the wild-type ShIR channel. In the KcsA channel, however, sensitivity to internal TEA is relatively low 5, and no TEA binding site in the inner vestibule, equivalent to position 441 located along the inner vestibule Shaker or to position 82 in the outer vestibule of KcsA, has been identified. For this reason, we do not use any special energy function for the intracellularly applied TEA.

As a TEA molecule carries one positive charge, one additional Cl⁻ for each molecule is placed in the internal reservoir to render the solutions in the reservoir electro-neutral. When one TEA molecule is placed in the reservoir, its concentration becomes 20mM. The concentration of K⁺, unless noted otherwise, is kept constant at 300mM. A TEA molecule enters the intracellular mouth of the channel quickly, attracted by the energy-well created by four glutamate residues, E₁₁₈. With low applied potentials, a TEA molecule that manages to enter the channel stays mainly in the energy-well in the intracellular gate, as shown in the dwell histogram shown in Figure 2A. With high applied potentials, however, the resident TEA molecule moves further inside the channel, and spends part of the time near the entrance of the selectivity filter (Figure 2BC). Because of its large size, a TEA molecule is unable to pass through the selectivity filter. In the presence of a TEA molecule anywhere along the ion-conducting pathway, the channel is effectively blocked, as no K⁺ ions can pass around this bulky molecule.

Display large version of this figure

Figure 2

Dwell histograms of internal TEA molecules in the channel. The inset shows the simulation assembly consisting of the KcsA channel at the center and two TEA molecules, K⁺ ions (pink), Cl⁻ ions (green) in or in the vicinity of the channel. The channel extends from −34.2Å to +34.2Å. The channel is divided into 100 thin sections and the average number of TEA molecules in each slice during a simulation period of 0.8μs is tabulated. The histograms are obtained with an applied potential of 0mV (A), +112mV (B), and +168mV (C).

The amount of blocking depends both on the applied potential and TEA concentration. In Figure 3A, the percentages of the currents relative to the control currents are plotted against applied potentials with the TEA concentrations of 20, 40, and 100mM. With the concentration of 20mM, we see that TEA blockade is weakly voltage-dependent when the applied potential is low, but becomes strongly voltage-dependent at higher applied potentials. With higher TEA concentrations, however, the currents relative to the control currents rapidly decline with an increase in the applied potential. The inhibitory constant for the internal TEA block we derive is 40mM at 130mV and 20mM at 170mV. The results described thus far are obtained from simulations, in which a TEA molecule is represented as a sphere with a radius of 3.4Å. The minimum and maximum radii of a TEA molecule are, respectively, 2.3 and 4.6Å, and its volume average is 3.4Å. To test the sensitivity of the results on the radius of sphere used, we repeat one set of simulations as shown in Fig. 3, representing a sphere with a radius of 2.3Å, instead of 3.4Å. Using a TEA concentration of 20mM, we determine the percentage of block, relative to the control, at various applied potentials. The results, shown in the uppermost curve in Figure 3A (open diamonds) may be compared with those obtained with a sphere of a radius of 3.4Å. There is a small, but consistent, decrease in the effectiveness of internal TEA blockade when the radius representing the molecule is reduced. This is because a K⁺ occasionally is able to pass around a TEA molecule blocking the intracellular gate of the channel.

Display large version of this figure

Figure 3

Reduction of the currents by internally applied TEA molecules at various applied potentials. (A) The current at each applied potential is expressed as a percentage of the current in the absence of TEA. The concentrations of TEA in the internal reservoir are 20mM (solid circles), 60mM (open circles), and 100mM (solid triangles). These three curves are obtained by representing a TEA molecule as a sphere with a radius of 3.4Å. The measurements are repeated, for the TEA concentration of 20mM using a sphere with a radius of 2.3Å (open diamonds). In this and all subsequent figures, error bars have a length of 2 SEM, and are not shown when they are smaller than the data points. (B) The current-TEA-concentration curve is obtained with the ionic concentration of 150mM and the applied potential of 224mV. Superimposed on the simulated data (solid circles) are the experimental measurements (open diamonds) obtained by Heginbotham et al. 5.

For KcsA, the experimentally determined value of the inhibitory constant k_i for internally applied TEA is 22mM 5. The applied potential and ionic concentration used to determine this constant are, respectively, 200mV and 100mM. To check how closely our results match the experimental observations, we construct a current-TEA-concentration for the ionic concentration curve, using the K⁺ concentration of 150mM and the applied potential of 224mV. The experimental measurements (open diamonds) are superimposed on the curve derived from BD simulations (solid circles) in Figure 3B.

A TEA molecule in the external vestibule of KcsA sees an attractive potential well created by four tyrosine residues at position 82, located just outside the selectivity filter. Molecular dynamics calculations reveal that this well extends ∼6Å along the z axis, as shown in Fig. 4. The one-dimensional free energy profile (dashed line) obtained by Crouzy et al. 9 is converted to one-dimensional potential energy (solid line). The resulting profile, once the free energy profile has been converted into the potential energy profile to be used in BD simulations, has a shape virtually identical to the free energy profile. The well starts at z=30Å (or z=15Å from the center of mass of the selectivity filter) and reaches its deepest point at z=24.44Å. The depth of the potential well at this position is 10 kT, ∼2.5 kT deeper than the one-dimensional free energy profile. We incorporate this potential energy into the BD algorithm, taking into account the short-range repulsive forces.

Display large version of this figure

Figure 4

Conversion of the one-dimensional free energy profile to a one-dimensional potential energy profile. The analytical function fitted to the experimentally determined free energy profile (solid line), as shown in Fig. 1, is used to derive an equivalent one-dimensional potential energy profile (dashed line) to be used in Brownian dynamics. The theoretical basis for this relationship is given in Appendix 1 . When a short-range force is added to this potential profile, it becomes equivalent to a three-dimensional free energy profile that can be incorporated into the Brownian dynamics algorithm.

Using this profile, with the well-depth of ∼10 kT at z=24.44Å, we construct the current-TEA-concentration profiles at three different values of applied potentials. For all potentials, the currents decrease exponentially as a function of the TEA concentration, as illustrated in Fig. 5. The inward currents are reduced to 50% of the control value at the TEA concentrations of 10, 17, and 24mM when the applied potentials are, respectively, −224, −168, and −112mV (inside negative with respect to outside). Thus, the magnitude of attenuation of the inward current by TEA depends on its concentration as well as the applied potential. The inhibitory constants k_i we obtain at all applied potentials are higher than the experimentally determined value, which is 3.2mM at 200mV 5.

Display large version of this figure

Figure 5

The current-TEA-concentration curves. The one-dimensional potential energy profile shown in Fig. 4, after adding a short-range force, is incorporated into the BD algorithm and the current flowing across the channel at different TEA concentrations is determined under the applied potentials of −224mV (solid triangles), −168mV (solid circles), and −112mV (open diamonds). The data points are fitted with the exponential function of the form: I=I₀ exp(−bc), where I₀ is the control current in pA and c is the TEA concentration in mM. The values of I₀ and b for the three curves are: 11.4 and 0.065 for −224mV; 7.2 and 0.043 for −168mV; and 4.36 and 0.037 for −112mV.

How deep does the energy-well have to be to replicate the experimentally determined inhibitory constant? We systematically change the depth of the potential energy profile by multiplying the original curve by 0.9, 1.1, 1.2, 1.3, and 1.4 and then construct the current-concentration curves with the applied potential of −168mV. The corresponding well-depths at z=24.44Å are, respectively, 9, 11, 12, 13, and 14 kT, as shown in Figure 6A. The family of current-concentration curves obtained with five different profiles of differing well-depths are displayed in Figure 6B. The curve obtained with the original profile (reproduced from Fig. 5) is included here for completeness. Each set of data points is fitted with an exponential curve. The half saturation values k_i obtained from the current-concentration curves are: 9.6, 5.4, 3.1, and 2.0mM when the original profile is multiplied by 1.1, 1.2, 1.3, and 1.4, respectively.

Display large version of this figure

Figure 6

The current-TEA-concentration curves with varying well-depths. (A) The depth of the three-dimensional free energy profile is systematically altered by multiplying the one-dimensional potential energy profile shown in Fig. 4 by 0.9, 1.0, 1.1, 1.2, 1.3, and 1.4. The depth of the original potential energy profile at z=24.44Å from the center of the mass of the selectivity filter is 10.0kT. (B) For each profile shown in panel A, the current-TEA-concentration curves are constructed. The free energy profiles used for simulations are indicated in the inset. The parameters b of the exponential function, I=a exp(−b(c)), used to fit the data points from the top to bottom curves, are 0.026, 0.047, 0.099, 0.194, 0.328, and 0.446.

We next explore how external TEA blockade is influenced by the applied potential. With a fixed TEA concentration of 3.2mM TEA, and with the potential profile with the well-depth of 13 kT, we determine the currents flowing across the channel at various applied potentials. The current-voltage relationships obtained with (open circles) and without TEA (solid circles) in the external reservoir are illustrated in Fig. 7. The currents under both conditions increase linearly with the applied potentials. The percentage of currents attenuated by 3.2mM TEA in the external vestibule, however, changes from 70% of the control current at the applied potential of −56mV to 30% at −256mV. Thus, for a fixed concentration, external TEA blockade becomes progressively more effective as the driving force is increased.

Display large version of this figure

Figure 7

Reduction of the currents at different applied potentials. The TEA concentration in the external reservoir and the depth of the well seen by a TEA molecule at z=24.44Å are kept constant at 3.2mM and 13kT, respectively. (A) The current increased linearly with an increasing applied potential in the absence (solid circles) and presence (open circles) of TEA. (B) The percentage of the current relative to the control current decreases steadily as the applied potential is increased.

The inhibitory constant k_i also depends on the concentration of KCl in the outer reservoir. We construct, in Fig. 8, the current-KCl-concentration curves with (open circles) and without TEA (solid circles) in the external reservoir. The TEA concentration and the depth of the potential profile are, respectively, 3.2mM and 13 kT, as used for the simulations illustrated in Fig. 7. The current under both conditions increases rapidly with an increasing ionic concentration and then saturates, leading to current-concentration relationships of the Michaelis-Menten form. The percentage of attenuation caused by TEA of a fixed concentration decreases as the ionic concentration is increased, as shown in Figure 8B. For example, the current in the presence of 3.2mM TEA is reduced to 31% of the control value when the ionic concentration in the reservoir is 75mM. The corresponding value at the KCl concentration of 900mM is 60%.

Display large version of this figure

Figure 8

Effectiveness of internal TEA at different K⁺ concentrations. The TEA concentration in the external reservoir and the well-depth seen by a TEA molecule at z=24.44Å are kept constant at 3.2mM and 13kT, respectively. (A) The currents increase both in the absence (solid circles) and presence (open circles) of TEA increase monotonically as the K⁺ concentration is increased from 75mM to 900mM. (B) The percentage of the current relative to the control current increases as the K⁺ concentration is increased.

The best agreement between the experimentally determined k_i value for KcsA 5 and simulation results is obtained when the depth of the potential profile seen by TEA is assumed to be 13kT. In Fig. 9, we plot the current-TEA-concentration curve obtained with the driving force of 196mV and the ionic concentration in both reservoirs of 150mM. Superimposed on the simulated data (solid circles) are the measurements obtained by Heginbotham et al. 5 (open diamonds). These experimental data points are obtained with the ionic concentration of 100mM and the applied potential of 200mV. To replicate the k_i value for the mutant Shaker T449Y 4, the depth of the free energy-well needs to be increased to ∼15kT.

Display large version of this figure

Figure 9

Comparison with the experimentally determined inhibitory constant. The currents derived from BD simulations (solid circle) at various TEA concentrations are compared with those determined experimentally (open diamonds) by Heginbotham et al. 5. The applied potential and the ionic concentration used for BD simulations are, respectively, 196mV and 150mV. The depth of the energy-well encountered by a TEA molecule is assumed to be 13kT.

The results from BD simulations show that externally applied TEA binds to the binding pocket located just exterior to the selectivity filter and becomes unbound according to a first order Markovian process. We tabulate the time the channel is not blocked by a TEA molecule (off-time) and the time the binding pocket is occupied by a TEA molecule (on-time). A short segment of the on-time and off-time sample path generated from BD simulations is illustrated in Fig. 10. When bound, a TEA molecule effectively blocks the passage of K⁺ ions. Conduction events occur only when no TEA is in the binding pocket, as indicated with downward marks in Fig. 10. With a 3.2mM TEA concentration, the channel is blocked and unblocked nearly 300 times in 100μs, which is usually one digitizing interval of patch-clamp recordings.

Display large version of this figure

Figure 10

Sequence of TEA binding and unbinding. A short segment of the TEA binding and unbinding process is illustrated. The sequence is governed by the first-order Markov process, and the conduction events take place only when a TEA molecule is not in the binding pocket, as indicated by downward marks.

To characterize the binding and unbinding processes, we carry out a set of simulations at eight different TEA concentrations under the applied potential of 168mV and calculate the mean off-time and on-time at each concentration. The depth of the well and ionic concentration used for these simulations are 12.0 kT and 300mM, respectively. The random sequence of off-times and on-times generated from the BD simulation can be statistically modeled as a first-order two-state continuous-time Markov process with a transition rate matrix denoted Q(c) as a function of the TEA concentration c. We statistically validate whether or not the binding or unbinding process obeys a first-order, continuous-time, Markov process. One way to demonstrate this would be to construct an interval histogram, which will be distributed as a single exponential function if the process is a first-order Markovian. That is, the random sequences of on-times or off-times denoted t for a given concentration c have an exponential probability distribution function of the form 1-exp(−t/τ_ON(c)) and 1-exp(−t/τ_OFF(c)), respectively. Here, we use the well-known Anderson-Darling statistical test 33 to determine whether the small number of random sequence of on- and off-intervals we obtained from BD simulations can be modeled by an exponential distribution. By using the Anderson-Darling test with a standard significance level of 5%, we find that the hypothesis that the sequence of off-times and on-times generated by BD simulation at various TEA concentrations is exponentially distributed, cannot be rejected. Thus, the random sequence of on-times (blocking of the channel by TEA) and off-times (unblocking of the channel) obtained from BD simulations forms a two-state continuous-time Markov chain holds.

The mean off-times τ_OFF(c) and on-times τ_ON(c) for different concentrations c in mM computed from the data is shown in Figure 11A. The figure shows that the time the molecule stays bound to the channel at different TEA concentrations, τ_ON, remains fairly constant (open circles), whereas the time channel remains unbound to the molecule, τ_OFF, decreases as the TEA concentration in the reservoir increases (solid circles).

Display large version of this figure

Figure 11

Concentration-dependent kinetics of TEA biding and unbinding. (A) The mean off-time τ_OFF(c) (solid circles) and the mean on-time τ_ON(c) (open circles) are plotted against the TEA concentration. The on- and off-durations are tabulated during a simulation period lasting 6.4μs. (B) The stationary probability P_OFF of a TEA molecule being not bound to the channel is plotted against the TEA concentration.

We now show that the elements of the transition rate matrix Q(c) behave as exponential functions of the concentration c. Thus, the binding and unbinding of TEA from BD simulations has a compact statistical representation. By relating the continuous-time transition rate matrix Q(c) to the TEA concentration c, we give an explicit characterization of the stationary distribution of the Markov chain and hence the average channel current versus concentration c. This allows us to characterize the effect of the blocker on the channel current and in particular give explicit trends as to how the average channel current from BD mimics the actual experimental channel current.

Let us denote the transition rate matrix of this two-state continuous-time Markov chain as

(9)

Our next task is to model how Q(c) evolves with concentration c. The stationary distribution of the continuous-time Markov chain, which gives the normalized times in which the channel is blocked (P_ON(c)) and not blocked (P_OFF(c)), is given by the solution of

(10)

(11)

A similar argument can be made for a random sequence of off-times and on-times as a function of the applied voltage. The depth of the well and the TEA concentration used for this series of analysis are, respectively, 13 kT and 3.2mM. In this case, off-time τ_OFF (solid circles) and on-time τ_ON (open circles) are plotted in Figure 12A against applied potentials. The average time TEA stays bound to the channel, once it enters the binding pocket, increases as the applied potential increases, whereas the average time the channel remains unbound from the molecule decreases with an increasing applied potential. Once again, all random sequences obtained from BD simulations are shown to be exponentially distributed using the Anderson-Darling test. In Figure 12B, the P_OFF(V) is fitted with the following two exponential functions: exp(−0.003V and −0.133), where V is in mV.

Display large version of this figure

Figure 12

Voltage-dependent kinetics of TEA binding and unbinding. (A) The mean off-time τ_OFF (V) (solid circles) and the mean on-time τ_ON (V) (open circles) are plotted against applied potentials. (B) The stationary probability P_OFF of a TEA molecule being not bound to the channel is plotted against applied potentials.

For this study, we approximated the shape of TEA as a sphere with a radius of either 2.3Å or 3.4Å. Ideally, the molecule should be represented as an oblate spheroid whose long and short axes correspond to 2.3 and 4.6Å. To incorporate spherically asymmetrical molecules as rigid bodies into the BD algorithm, we would need to follow, in addition to their translational steps, the rotational motions using unbiased operators for finite rotations 34. We should also need to determine the diffusion coefficient for rotational motions and may have to incorporate a hydrodynamic interaction term into the diffusion tensor 35. For this preliminary study, however, we use a simplified, spherical TEA model, deferring the more realistic representation of the molecule to a future project.

The stabilization of TEA by the four tyrosine residues in KcsA is likely to arise either from the hydration forces 9,10 or from the attractive interaction between a TEA molecule and the quadrupolar moment of the aromatic ring 3,11. Because in Brownian dynamics water molecules are not explicitly simulated, the interaction originating from the hydration force is not taken into account. Similarly, forces stemming from electrostatic interactions between a charged particle and quadrupolar moments are not calculated in our algorithm. We therefore implicitly account for these forces by incorporating in the BD algorithm the free energy profile obtained from molecular dynamics, so that each TEA molecule in the reservoir sees, in addition to all other attractive and repulsive forces it encounters, the attractive force generated by the tyrosine residues near the entrance of the selectivity filter. The so-called potential of mean force derived from umbrella sampling calculations using the weighted histogram analysis method 36 is in fact a one-dimensional free energy profile. For use in BD, it should first be converted into a one-dimensional potential energy profile and then into a three-dimensional free energy profile. The interrelationships between these profiles are briefly summarized in Appendix 1 .

By grafting into BD the additional force encountered by TEA in the simulation assembly, we uncover several subtle features of its mechanism of blockade. First, we determine the precise depth of the free energy profile external TEA has to encounter to reproduce the experimentally determined inhibitory constant (Figure 6 and Figure 9). We show that the depth of the energy-well has to be somewhat deeper than that calculated from molecular dynamics 9 to account for the experimental measurements. Secondly, the affinity to external TEA block also depends on the applied potential (Fig. 7) and, to a lesser extent, the ionic concentration in the outer reservoirs (Fig. 8). At a fixed applied potential, the current decreases exponentially with an increasing TEA concentration. Likewise, with a fixed TEA concentration, the percentage of inhibition relative to the control current increases as the applied potential increases. Thus, an inhibitory constant k_i derived at a fixed potential with one ionic concentration may not hold true for different experimental conditions, as noted previously by Eisenberg 37. Thirdly, our analysis of the sequences of blocking events reveals that binding and unbinding of TEA to the binding pocket is a first-order Markovian process. The transition rate for TEA binding increases steeply with the concentration, whereas that for unbinding remains invariant of the concentration (Fig. 11). The transition rate for binding and unbinding depends similarly on the applied potential (Fig. 12).

The depth of the well from the free energy profile calculated from molecular dynamics is 4.5 kcal/mol or 7.5 kT 9. The binding energy computed by Luzhkov and Åqvist 10 using an automated docking approach is slightly lower, ranging from 3.3 to 3.9kcal/mol. From the one-dimensional free energy profile, we can directly estimate the probability of a TEA molecule found in the binding site, P_bind, and this in turn leads to an estimate for k_i (see Appendix 2 ). The value of the inhibitory constant with no applied potential estimated from the free energy profile of Crouzy et al. 9 is 23mM. When the applied potential is increased from 0 mVA to 200mV, the value of k_i shifts by 2.1mM. Although the free energy profile obtained from molecular dynamics is unable to account for other subtle features of the dynamics of TEA blockade, it nevertheless provides an order-of-magnitude estimate of the inhibitory constant. More importantly, such a profile gives valuable information about the general shape of the three-dimensional free energy profile that needs to be incorporated into the Brownian dynamics algorithm. It is possible that the depth of the free energy profile calculated from molecular dynamics may become deeper once all four tyrosine residues are repositioned such that a TEA molecule in the binding site would interact with the faces of the aromatic side chains. Such an en face orientation is needed for the favorable cation-π electron interaction 11.

Yellen et al. 2 noted that internal TEA blocks a mutant Shaker channel, ShIR, at a submillimolar concentration. When threonine residues at position 441 are replaced with serines, the sensitivity to internally applied TEA is reduced by an order of magnitude. It thus appears that the Shaker K⁺ channel possesses an internal TEA binding pocket similar to that found in the external vestibule. In contrast, the inhibitory constant for internally applied TEA for the KcsA K⁺ channel is 22mM 5. This value is close to the results of our simulation (Fig. 3). The suppression of the currents by internal TEA in our simulations is caused by the occlusion of the ion conducting pathway. The blocking molecule penetrates deeper inside of the pore as the applied potential is increased, thus making the effectiveness of blockade highly voltage-dependent (Fig. 2).

By making several simplifying assumptions, we have been able to deduce many of the salient features of internal and external TEA blockade in the KcsA K⁺ channel. Here we have exploited the capability of BD to measure the current flowing across a biological ion channel to study how internally or externally applied TEA interacts with the permeation dynamics. By carrying out BD simulations in the KcsA potassium channels with TEA ions in KCl solutions, we show how the inward and outward currents change with an increasing concentration of TEA in the simulation assembly. From a series of such current-concentration curves constructed from different energy profiles with varying well-depths, we are able to infer the depth of the energy-well a TEA molecule needs to see to replicate the experimentally determined inhibitory constant for the TEA block. We also demonstrate that the attenuation of currents by TEA is dependent on other factors, such as the applied potential and ionic concentrations of potassium ions in the reservoirs. Moreover, by tabulating how long a TEA molecule stays bound in the binding site, we reveal the dynamics of binding and unbinding processes taking place in a fast timescale. The technique we introduce here can be fruitfully utilized in examining other K⁺ channel models whose atomic coordinates are not yet determined. Any homology model of a K⁺ channel constructed must account for its current-voltage-concentration profiles as well as the effects of quaternary ammonium ions on channel currents.

Appendix 1

one-dimensional and three-dimensional free energy profiles

The aim below is to determine the relationship between the three-dimensional potential model that can be used in BD and the one-dimensional free energy profile obtained by using molecular dynamics calculations 9. In BD, we model the three-dimensional potential as

(12)

where U^PMF(z) is an arbitrary one-dimensional potential, and U^SR(x, y, z) models the short-range forces due to the channel walls. We wish to determine the potential U^PMF(z) so that the one-dimensional energy profile corresponding to Eq. (12) is equal to the one-dimensional free energy profile determined by molecular dynamics. As the free energy profile already includes an entropic component, we need to calculate the equivalent entropic component in BD. To do this, we must find the probability that the TEA resides within each cross-sectional area of the KcsA channel.

We start with the one-dimensional free energy profile obtained by using molecular dynamics. In one dimension, and at equilibrium, the probability of finding an ion per unit length is

(13)

where G_1d(z)=G^PMF(z) is the free energy in one dimension taken from Crouzy et al. 9, and β=1/kT. Here, k_1d is a normalization constant, and L₀ is a reference length.

Even though profiles generated by WHAM analysis of molecular dynamics simulations are described as potentials of mean force, the nature of the analysis means that they reflect the probability of finding an ion at positions along the reaction coordinate, rather than the average (systematic) force along the reaction coordinate. In other words, they are free energy profiles, not potential energy profiles, and include an entropic component which depends on the effective cross-sectional area,

(14)

where U^PMF(z) is the potential energy, T is the absolute temperature, and −TS^PMF is the entropic component of the free energy.

In the case of hard boundaries, this entropic component is given by

(15)

where A^PMF(z) is the cross-sectional area at each z position and A₀=(L₀)² is a reference area. A larger area gives a larger entropy which leads to a smaller free energy and a greater probability (per unit length) of occupancy by an ion.

Now consider the three-dimensional potential model that can be used in BD. The probability of finding the TEA molecule per unit volume for any given cross section is

(16)

where G_3d=U^BD(x, y, z) is the free energy in three dimensions and k_3d is a normalization constant. Again, this is only valid at equilibrium, but the process of TEA binding and unbinding occurs at equilibrium.

The one-dimensional and three-dimensional probabilities can be related using

(17)

By substituting from Eqs. (13), this gives

(18)

After applying the natural log, Eq. (18) becomes

(19)

Since G_3d=U^BD(x, y, z) from Eq. (12), Eq. (19) becomes

(20)

After rearranging terms,

(21)

The integral is an effective cross-sectional area, taking into account the short-range repulsive forces, rather than a hard boundary.

We assume that the channel has cylindrical symmetry such that

(22)

This assumption is justified since U^SR, the short-range forces from the BD force field, have cylindrical symmetry. Because G_1d(z) is the result from Crouzy et al. 9, Eq. (22) represents the required correction needed to solve for the U^PMF(z) for use in BD simulations.

Eq. (22) can be rewritten as

(23)

where

(24)

is the dimensionless free-energy offset, and

(25)

is the effective cross-sectional area, taking into account the short-range forces.

To make use of Eq. (23) in a BD simulation, we need to set a value for k_G so that the potential derived from the profile derived from molecular dynamics will match up with the potential from macroscopic electrostatics at the limit of the potential of mean force. We assume that at z_max, the extracellular limit of the free energy profile,

(26)

so that, from Eq. (23),

(27)

Appendix 2

estimating the inhibitory constant from pmf

It is possible to compute P_bind, the probability of a TEA being found in the binding site, directly from the one-dimensional free energy profile and the concentration of TEA in the reservoir. This in turn leads to an estimate for k_i, the concentration at half inhibition, by setting P_bind=1/2.

We can estimate

(28)

(29)

Since G_1d(z)≡G^PMF(z), all we need in addition to the free energy profile to estimate P_bind from this formula is an estimate of k_1d.

From Eq. (23), we obtain

(30)

But U^PMF(z_max)−G_1d(z_max)=0 from Eq. (26), so

(31)

We take k_3d/V₀=k_RES, where V₀=(L₀)³, to be the TEA concentration in the reservoir in appropriate units (number density). Then from Eqs. (29)

(32)

If P_bind=1/2, then k_i=k_3d/V₀ and

(33)