INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Gramicidin A (GA) is an antibiotic polypeptide that consists of 15 amino acid residues. Its -helical, head-to-head dimer in membranes forms a single-file ion channel (Urry, 1971) that selectively conducts monovalent cations, binds divalent cations, and rejects all anions (for reviews see Andersen and Koeppe, 1992; Busath, 1993; Koeppe and Andersen, 1996; Wallace, 1998). Recent high-resolution structures of the GA channel have been determined from solution NMR (Arseniev et al., 1985) and solid-state NMR studies (Ketchem et al., 1993, 1997; Separovic et al., 1999). In the absence of structural information for biological ion channels, the GA channel has been the main focus of theoretical investigations for a long time (Partenskii and Jordan, 1992; Roux and Karplus, 1994). The recent determination of the crystal structure of the KcsA potassium channel (Doyle et al., 1998) has now shifted the attention to biological ion channels. Nevertheless, with its simple and well-defined structure and ample functional data, the GA channel should continue to play a prominent role as a test bed for theories of ion permeation.

Modeling of the GA channel has evolved from simple electrostatic calculations with rigid dielectric boundaries (Levitt, 1978; Jordan, 1982; Monoi, 1991) to sophisticated all-atom molecular dynamics (MD) simulations with GA embedded in a lipid bilayer and solvated with water (Woolf and Roux, 1994, 1997; Chiu et al., 1999). These developments are reviewed in several articles to which we refer for a complete list of references (Pullman, 1987; Partenskii and Jordan, 1992; Roux and Karplus, 1994). Most of the potential or free-energy profiles of ions obtained in these studies are in qualitative agreement with the observed binding sites at each end of the channel and a central barrier in between them. However, the calculated energy barriers for the translocation of ions are in the range of 20-30 kT (Jordan, 1990; Roux and Karplus, 1993), which are too high to predict the experimental conductances. Indeed, these profiles have been tested by McGill and Schumaker (1996) by inserting them into a diffusion theory of permeation. They found that the magnitudes of the profiles had to be reduced substantially to reproduce the observed currents. For the microscopic models, the problem with the profiles most likely lies in the force fields used in the MD simulations. Recent high-resolution NMR studies of cation transport in the GA channel (Tian et al., 1996; Tian and Cross, 1999) have shed further light on this problem by demonstrating that the GA peptide remains rigid upon cation binding and the ion is solvated by just two carbonyl oxygens. In contrast, early microscopic models typically predicted a rather flexible backbone with the carbonyl oxygens swinging by 20-40° so that four carbonyls and two water molecules provide adequate solvation for the cation in the channel. More recent studies suggest that a permeating sodium ion moves off axis to be solvated by the carbonyl oxygens so that the channel may be less flexible. Nevertheless, the carbonyl oxygens are still observed to rotate in substantial amounts in the presence of a sodium ion (Woolf and Roux, 1997). In the rigid scenario, the missing solvation energy is presumably provided by the water molecules in the channel having a more ordered structure than in the flexible models. In this regard, inclusion of the polarization effects in the standard MD force fields would be very desirable.

This development is a mixed blessing for the lower-resolution permeation theories such as the Poisson-Nernst-Planck (PNP) equations (Eisenberg, 1996, 1999) and Brownian dynamics (BD) simulations (Cooper et al., 1985; Chung et al., 1998, 1999). On one hand, the NMR experiments appear to justify the assumption of a rigid channel structure used in these theories, which is often one of the main criticisms of them. On the other hand, the long-range order of the water molecules in the channel imposed by this rigid structure creates even more serious problems for using continuum electrostatics in PNP and BD approaches. In a bulk solution, the electric field of an ion and hence the polarization of the surrounding water molecules drops as 1/r². Because of the focusing of the electric field by the dielectric boundary, the reduction in the field of an ion in a channel is not as severe as in bulk. Nevertheless, continuum electrostatics predicts substantial reduction in polarization of channel waters with distance from an ion, in disagreement with the MD predictions that the water dipoles are well aligned along the channel axis. Notwithstanding such criticism from the microscopic quarters (Partenskii and Jordan, 1992; Roux and Karplus, 1994; Roux et al., 2000), the PNP theory has recently been applied to the GA channel, giving an apparently successful description of the current-voltage relations (Kurnikova et al., 1999; Cardenas et al., 2000; Hollerbach et al., 2000). In view of the recent queries on the validity of the PNP theory in narrow channels (Corry et al., 2000a, b), this success is even more remarkable and warrants closer scrutiny. Thus, the main thrust of this work is to perform continuum electrostatic calculations using the high-resolution structure of the GA channel (Ketchem et al., 1997) and use these results in BD simulations to see whether these methods can provide a consistent description of the available experimental data on the GA channel.

There are several outstanding issues in the GA channel that we will attempt to address in this study. The most pressing question is the value of the effective dielectric constant in the channel, _c, and whether such a uniform value exists at all. While semi-microscopic calculations suggest _c = 3-5 (Partenskii et al., 1994), the bulk value of 80 is used in the PNP calculations quoted above. To resolve this question, the whole range of _c values will be considered when solving Poisson's equation. A second issue is the origin of cation selectivity of the GA channel. It has been argued that because the GA peptide is charge-neutral, there is no obvious mechanism to explain its valence selectivity. Various mechanisms based on intricate ion-peptide-water interactions have been proposed to explain this property (Sung and Jordan, 1987; Dorman et al., 1996; Roux, 1996). However, in all these semi-microscopic and microscopic calculations the GA peptide has a flexible structure, and it is not clear how many of these results would carry through if it were rigid. In addition, the PNP calculations of potential and concentration profiles inside the channel suggest that the valence selectivity can be understood simply in terms of the charge distribution in the peptide. It is important to ascertain whether the electrostatic interaction between an ion and GA peptide could indeed discriminate between cations and anions. Such a simple basis would obviate the need to look for more complicated explanations for valence selectivity.

If one cannot trust the free energy profiles obtained from MD simulations, and continuum electrostatics fails to give reasonable potential energy profiles, what other avenues are available for progress? Solving the inverse problem, that is, going from experimental data to potential may provide a more rewarding alternative to direct studies of the GA channel under the present circumstances. Using insights from experiments and simulations, one can construct a potential profile that provides a satisfactory description of physiological data when fed into BD simulations. Such a study was undertaken earlier by relating potential profiles to sodium currents using electrodiffusion equations and one-dimensional BD (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989). This inverse technique is extended here to potassium permeation using three-dimensional BD. Although one-dimensional BD should do a reasonable job of representing the single-file motion in the narrow GA channel, the extension to three dimensions has the advantage that the entry and exit of ions to the channel can be explicitly modeled. It also allows ions to move off-axis to interact more intimately with the carbonyl oxygens, as seen in MD simulation of sodium ions. Compared to the continuum theory, BD accounts for individual ion motions and self energies in a direct manner in which the electrodiffusion equations cannot.

Besides providing a better understanding and a united explanation for various channel observables, potential energy profiles found with this "inverse" technique would also be useful in related studies of the GA channel and its analogs. For example, it could be used in constraining the free-energy calculations in MD simulations when searching for more appropriate force fields in a channel environment. Another area rich in applications is the structure-function relationships in modified GA channels (Koeppe and Andersen, 1996). Changes in the GA structure have ranged from mutations in the amino acid sequence to fluorination of specific residues that modulate dipole strengths (Oiki et al., 1994; Andersen et al., 1998; Busath et al., 1998; Cotten et al., 1999; Borisenko et al., 2000; Anderson et al., 2001). Especially in the latter case, where changes in dipole strength can be easily incorporated in the potential profile for the native structure, one could predict the expected functional changes with a minimal computational effort. Here we construct a potential profile for potassium ions that could be used for such purposes in modified GA channels.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Gramicidin A structure

A high-resolution structure of the GA channel has recently been determined from solid-state NMR spectroscopy via refinements of the initial structure against both the experimental constraints and the CHARMM global energy (Ketchem et al., 1997). Here we use the atomic coordinates obtained by Ketchem et al. (1997) and stored in the Protein Data Bank with the accession code 1MAG. This structure is slightly different from that of Arseniev et al. (1985), which has been used in most of the GA model studies in the past (see below for a comparison of the two structures). For consistency with the high-resolution structure, the partial charges on the atoms are taken from the all-atoms PARAM22 (Mackerell et al., 1992) version of the CHARMM force field. The effect of using partial charges from another force field (AMBER, Cornell et al., 1995) will be discussed in the Results section.

To define the surface of the channel (and hence the dielectric boundary), the radii of atoms in the GA peptide are required. We use the minimum set derived by Li and Nussinov (1998, Table 6). Following their recommendation, the Coulombic radius is used for the polar N and O atoms and the van der Waals radius for the rest of the nonpolar atoms. (Had we used the van der Waals radius for all the atoms, the channel radius would be ~10% smaller, leading to a slightly higher self-energy for ions.) Slices of channel profiles generated by this data set are shown in Fig. 1, A and B.

View larger version (11K): [in this window] [in a new window]   FIGURE 1   (A) and (B) Two slices of the channel boundary down the central axis at different azimuthal angles . The regions of different dielectric constant are noted in (B). The values shift from the channel value (c) to the bulk value (b) smoothly over a 5 Å distance centered about the dotted lines. (C) A transverse section of the axially symmetrized channel shape.

The GA dimer is embedded in a neutral membrane of length 33 Å, modeled as a uniform dielectric medium (with dielectric constant equal to that of the GA peptide) without any charges or dipoles. This length and model matches that of relatively neutral glycerylmonoolein (GMO) membrane. It is also close to the thickness of diphytanoylphosphatidylcholine (DPPC) (35 Å), so that results of BD simulations can be compared with experiments involving both types of membranes, although the dipoles in the DPPC membrane may affect the current. Because GMO is wider than the GA dimer, hydrophobic matching is used in smoothly joining the membrane to the GA headgroups (Harroun et al., 1999) as indicated in Fig. 1. Finally, to complete the simulation system, cylindrical reservoirs of 30 Å radius and length are connected to each end of the channel. When filled with an electrolyte solution, these reservoirs provide an adequate representation of the intra- and extracellular baths.

As described below, the calculation of potentials is too slow to be carried out at each time step during BD simulations, and is too cumbersome to be precalculated and stored when using a complete three-dimensional shape. For this reason, in BD simulations, an axially symmetric form of the channel boundary is used for convenience (Fig. 1 C). This shape is generated by averaging the radial distance around the circumference and modifying it further until the calculated potential profile roughly matches the one obtained from the original shape shown in Fig. 1, A and B. Provided these potential profiles are very similar, the use of a symmetric shape should not modify our results.

Solution of Poisson's equation

In previous BD studies of channels we have used a boundary element method to solve Poisson's equation (Levitt, 1978; Hoyles et al., 1998). This method works fine when the channel and reservoir waters have the same dielectric constant (_c = _w = 80), but becomes problematic when they differ. For small reductions in _c values (e.g., _c = 60), one can use the smaller value everywhere and incorporate the neglected Born energy as a barrier at either entrance of the channel without much loss in accuracy (Chung et al., 1999). However, for larger reductions in _c values one needs to use a finite difference method that allows such a change to be included explicitly. For this purpose, we have refined the Poisson-Boltzmann solver developed earlier (Moy et al., 2000) by smoothing over the charge distribution and dielectric boundary (Davis et al., 1991). This code (with zero electrolyte concentration) is used in calculating the potential profiles for various ions types in the GA channel.

We have subjected this finite difference code to several tests to check its convergence properties and accuracy. In these tests _c = 80 is used for ease of comparison with the boundary element method. Fig. 2 A shows the grid size dependence of the potential energy profile of a monovalent cation as it is moved along the central axis. As the grid size in the finite difference method is reduced from 0.8 Å to 0.4 Å, the results are seen to converge rapidly. In Fig. 2 B we compare the potential energy profile obtained from the finite difference solutions (dashed line) with that obtained from the boundary element method (solid line). Here the axially symmetric channel boundary is used in both methods. A general agreement is obtained between the two methods (differences are much smaller compared to the thermal fluctuations), which establishes the validity of the finite difference solutions.

View larger version (15K): [in this window] [in a new window]   FIGURE 2   (A) Convergence of the potential energy profile with the grid size used in the finite difference solution of Poisson's equation. (B) Comparison of the energy profile obtained from the finite difference solution (dashed line) with that obtained from the boundary element method (solid line) for the axially symmetric boundary; c = 80 and p = 2 are used in both cases.

Brownian dynamics

A review of Brownian dynamics simulations in ion channels is given by Cooper et al. (1985). The one-dimensional channel models used in the early studies have recently been extended to realistic three-dimensional channel geometries. In particular, BD simulations have been applied in studies of model toroidal (Li et al., 1998), acetylcholine receptor (Chung et al., 1998), KcsA potassium (Chung et al., 1999, 2002), porin (Schirmer and Phale, 1999; Im et al., 2000), and L-type calcium (Corry et al., 2001) channels. We give a brief description of the method here and refer to the earlier work for details.

In BD, the motion of individual ions is simulated using the Langevin equation:

(1)

where m_i, q_i, and v_i are the mass, charge, and velocity of the ith ion. In Eq. 1, the effect of the surrounding water molecules is represented by an average frictional force with a friction coefficient m_ii, and a stochastic force F_R arising from random collisions. The last two terms in Eq. 1 are, respectively, the electric and short-range forces acting on the ion. The total electric field at the position of the ion is determined from solution of Poisson's equation, and includes all possible sources due to other ions, fixed and induced surface charges at the channel boundary, and the applied membrane potential. Because solving Poisson's equation at each time step is computationally prohibitive, we store precalculated values of the electric field and potential due to one- and two-ion configurations in a system of lookup tables, and interpolate values from these during simulations (Chung et al., 1999). The tables require a large amount of computer memory, and storage of three- and six-dimensional tables for a general boundary creates problems, even on a supercomputer. To facilitate BD simulations we generate an axially symmetric boundary that reproduces the original potential profile, and thereby reduce the dimensions of the tables by one. At short ranges, the Coulomb interaction between two ions is modified by adding a damped oscillating potential that replicates effects of the short-range repulsive and hydration forces (Corry et al., 2001).

The Langevin equation (1) is solved at discrete time steps following an algorithm devised by van Gunsteren and Berendsen (1982). To simulate the short-range forces more accurately we use a multiple time-step algorithm in our BD code. A shorter time step of 2 fs is used across the channel where short range ion-ion and ion-protein forces have the most impact on ion trajectories; elsewhere, a longer time step of 100 fs is used. If an ion is inside the short time step region at the beginning of a 100-fs period, then that ion is simulated by 50 short steps while the other ions in the long-time regions are frozen to maintain the synchronicity. Simulations under various conditions, each lasting for one million time steps (0.1 µs), are repeated numerous times. Initially, a fixed number of ions are assigned random positions in the reservoirs with velocities assigned according to the Maxwellian distribution. The current is determined from the number of ions traversing the channel during the simulation period. To maintain the specified concentrations in the reservoirs, a stochastic boundary is applied: when an ion crosses the channel, say from left to right, an ion of the same species is transplanted from the right reservoir to the left. More details of the system boundaries and applied potential can be found in Corry et al. (2002).

Due to the narrowness of the gramicidin channel, ions and water molecules must move in single file. This presents a further constraint on the motion of ions in the channel not encountered in previous BD simulations. If there are two ions in the channel, then they will be separated by an integer number of "trapped" water molecules and the distance between the ions must remain relatively constant. The intervening water molecules will prevent the ions getting closer, and gaining a larger separation would create a vacuum in the channel. We ensured that this condition is held in our simulations by subjecting both ions to the following additional force whenever a second ion entered the channel:
Here r is the distance between the ions and r₀ is a reference separation equal to the nearest integer number of water molecule diameters when the second ion first appears within the narrow section of the channel. This allows the ions to move freely back and forth while constraining the separation between them to remain roughly constant while they are both in the pore. The choice of the water diameter in determining r₀ is not of critical importance, and we simply use the canonical value of 2.8 Å.

The BD program is written in FORTRAN, vectorized, and executed on a supercomputer (Fujitsu VPP-300 or Compaq AlphaServer SC). The time to complete the simulations depends on how often ions enter the short time-step regions. With 48 ions in the system, the CPU time needed to complete a simulation period of 1.0 µs (10 million time steps) is ~30 h. A temperature of 298 K is assumed throughout and a list of the other parameters used in the BD simulations is given in Table 1 (note that the diffusion coefficient is related to the friction coefficient m in Eq. 1 by the Einstein relation, D = m/kT).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Electrostatic calculations

Potential energy profiles of single ions reveal the binding sites and barriers in the channel and thus provide clues to the permeability characteristics of a proposed model. The absence of a well at an observed binding site or presence of a large barrier (which would prevent conduction) would be sufficient grounds to reject a model without carrying out expensive simulation work. The profiles in this section are obtained from a finite difference solution of Poisson's equation using the nonsymmetric channel boundary (Fig. 1, A and B). In Fig. 3 A we show how the potential energy profile of a monovalent cation in the GA channel changes as the effective dielectric constant of channel waters is reduced from _c = 80 to 5. Here the dielectric constant of the GA peptide is set to _p = 2, representing its electronic polarizability. The very low values of _c ~ 5 suggested by microscopic calculations (Partenskii et al., 1994), are seen to lead to huge barriers (~55 kT). This will prevent conduction of ions at any realistic applied potential, as also noted by Partenskii et al. (1994). When a larger polarizability is assumed for the GA peptide (_p = 5), the barriers are reduced by a factor of 2-3 but still remain too large to permit conductance for low _c values (Fig. 3 B).

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Dependence of potential energy profiles on the effective dielectric constant of water in the channel (c). The dielectric constant of the protein is p = 2 in (A), and p = 5 in (B).

At the other extreme of _c = 80 used in the PNP calculations (Kurnikova et al., 1999; Cardenas et al., 2000; Hollerbach et al., 2000), the energy profile is almost flat, as shown in the lowermost profile in Fig. 3 A. (For the detailed features of this profile, see Fig. 2 A). At first glance, this profile with wells near the entrances and a central barrier appears quite reasonable. However, at only 1.5 kT, the wells are not deep enough to provide binding sites, nor is the barrier enough of an impediment to lead to the saturation of current. These points will become obvious when we present concentration profiles and current-concentration curves obtained from BD simulations below. In contrast, the potential profiles obtained in PNP calculations exhibit a deep potential well across the entire length of the channel. The discrepancy has been shown to arise from the neglect of ion self-energies in PNP (Corry et al., 2000b). To make this point more explicit, we decompose the profile obtained with _c = 80 and _p = 2 into the self-energy part due to the reaction field from the dielectric boundary (calculated by setting the partial charges in the peptide to zero) and the ion-peptide interaction part due to the partial charges, as shown in Fig. 4 A. It is seen that the flat profile follows from the near cancellation of these two large components, which have opposite signs. A deep potential well would result if the self-energy term were ignored.

View larger version (16K): [in this window] [in a new window]   FIGURE 4   (A) Individual contributions to the energy profile for c = 80 and p = 2 (middle) from the self-energy (top) and ion-peptide interaction terms (bottom). (B) Potential energy profiles for a divalent cation (solid line) and a monovalent anion (dashed line) for the above choice of  values. The monovalent cation result (dotted line) is included for reference.

Values of _c between the two extremes discussed above do not provide any improvement, either. While the barrier height increases with decreasing _c (see Fig. 3) the wells nearly disappear, thus disagreeing with the observed binding sites. This point will again be made clearer with the BD simulations below. In the following continuum electrostatic calculations we will use the common _c = 80 and _p = 2 values, as we have seen that variations from these only lead to a worse agreement with the experimental data.

We next consider the potential energy profiles for divalent cations and monovalent anions (Fig. 4 B). Experimentally, the former bind and block the GA channel and the latter are rejected. From a continuum electrostatics viewpoint, comparing the energy profiles provides a stringent test of the model because they are not independent quantities, but follow directly from the monovalent ion results in Fig. 4 A. In fact, if we denote the valence of the ion by z, the self-energy term of a monovalent ion by U_s, and the ion-peptide term by U_p, the profile for a divalent cation is given by U₂ = 4U_s + 2U_p and the one for monovalent anion by U₁ = U_s  U_p. Because U_p  U_s, one expects U₂  U₁  2U_s. This explains the similarity of the profiles for divalent cations and monovalent anions in Fig. 4 B. The absence of any potential wells in the divalent cation profile provides the most direct and clear evidence for the failure of continuum electrostatics in the GA channel. When such a large barrier is used in BD simulations it results in negligible ion concentrations inside the channel, in contradiction with the experimentally observed binding of divalent cations (see below). The large barrier for anions, however, provides an obvious mechanism for valence selectivity of the GA channel. Just like the divalent cations, anions would not enter the channel in the presence of a 25 kT barrier.

One may question the reliability of such an inference on valence selectivity from a continuum electrostatics calculation; after all, we have just condemned its use in the GA channel. To answer this query, we need to consider the results in Fig. 4 in more detail. From microscopic calculations we know that both the self-energy and ion-peptide interactions are largely underestimated in electrostatic calculations with _c = 80. The large barriers observed at lower _c values would be canceled by the ion-dipole interactions due to the ordered channel waters, which is ignored in the continuum calculations. Clearly, these interactions are even more important for divalent cations to provide the permanent binding sites at the channel entrance. Switching the sign of a monovalent cation, however, basically changes the sign of the ion-peptide interaction, leaving the other terms more or less intact. If the electrostatic results in Fig. 4 are underestimates, as indicated by microscopic calculations, then the barrier for valence selectivity can only go up in more realistic calculations. Thus the proposed valence selectivity mechanism based solely on the distribution of charges in the GA channel appears to be a robust result.

It is important to check the sensitivity of the above results against variations in the GA structure. Two significant parameters in this respect are the coordinates and magnitude of the partial charges of the peptide atoms. In Fig. 5 A we compare the profile obtained using the structure of Arseniev et al. (1985) with that of Ketchem et al. (1997) shown in Fig. 2 A. The new profile has no wells and a 4 kT barrier, and could only lead to worse agreement with experimental data than the old one. Not surprisingly, the energy-minimized Ketchem structure leads to a lower profile than the Arseniev one. In Fig. 5 B we show the effect of changing the partial charges from the CHARMM set (solid line) to AMBER (dashed line). The two sets have similar charges for the backbone atoms, but the CHARMM charges are a factor of 2 to 3 larger for the side-chain atoms. This explains the relative increase in the AMBER profile at the entrances (less negative charges) and decrease at the middle (less positive charges). Although the change in the profile is not significant, the disappearance of the central barrier can be viewed as a backward step. It is also worth noting that the failure to reproduce divalent binding sites is likely to be less dependent on the precise atomic structure than the monovalent profiles. As the divalent profile is dominated by the self-energy, large variations in the position of the boundary would be required to obtain binding sites.

View larger version (15K): [in this window] [in a new window]   FIGURE 5   (A) Comparison of the potential energy profiles obtained using the structure of Ketchem et al., 1997 and Arseniev et al., 1985. CHARMM partial charges are used in both cases. (B) Effect of using different sets of partial charges on the potential energy profiles; CHARMM (solid line), AMBER (dashed line). The Ketchem et al. (1997) structure is used in both cases.

Brownian dynamics simulations

The potential energy profiles presented above explain in qualitative terms why continuum electrostatics fails in the GA channel. To make contact with experimental data and thus demonstrate this inadequacy more quantitatively, we next perform Brownian dynamics simulations. For reasons stated in the Materials and Methods section, the axially symmetric boundary is used in solving Poisson's equation here.

There have been many experimental I-V measurements in the GA channel, and the success of the PNP equations in matching these has been used as an argument for its applicability in narrow channels (Kurnikova et al., 1999; Cardenas et al., 2000; Hollerbach et al., 2000). Here we wish to point out that, as long as the diffusion coefficient of ions in the channel is a free parameter, fitting linear I-V curves poses no problem for any channel model. As an example, the I-V curve obtained from BD simulations of a 500 mM KCl solution with _c = 80 (filled circles) is compared to experimental data (O. S. Andersen, private communication) (open squares in Fig. 6). The bulk values of the diffusion coefficients are used for both ionic species in the reservoirs, but the diffusion coefficient of K⁺ ions inside the channel, D, is reduced by 80% to fit the data. Note that such a reduction is not required for the Cl ions because they do not enter the channel. It is worthwhile to emphasize that unlike PNP equations, the current in the channel is not linearly proportional to the diffusion coefficient in BD simulations. This shape of the current versus diffusion coefficient depends on the depth of the energy well and height of the energy barrier (see Fig. 11). For the energy profile used to obtain Fig. 6, the current is reduced by less than a factor of 2 when D is suppressed to 10% of its bulk value. A further reduction causes a steep decrease in current. It is worth noting that ions never pass each other in the channel during simulations, the single file nature of GA permeation is reproduced.

View larger version (11K): [in this window] [in a new window]   FIGURE 6   Current-voltage relationships obtained from BD simulations of a 500 mM KCl solution (filled circles fitted with a line). Error bars have length one standard error. Experimental data from O. S. Andersen (private communication) are shown by open squares.

The results plotted in Fig. 6 show that this model has no difficulty in fitting the I-V data of the GA channel. However, we next show that this continuum electrostatic model of the GA channel fails completely when its predictions are compared with the observed binding sites and conductance-concentration curves. In Fig. 7 we show the concentration profiles of K⁺ ions in the channel without an applied potential (A), and with a 200 mV applied potential (B) obtained from BD simulations using a 500 mM KCl solution as above. The potassium concentration in the channel displays little hint of the expected binding sites at around ±11 Å (Tian and Cross, 1999). To be consistent with these data we should see two large concentration peaks separated by an obvious depression. In our BD data the central barrier causes a dip in the concentration profile in the absence of a driving force (A), but as soon as a potential difference is applied, this dip disappears (B), as ions can easily climb it. Thus, both the wells and the barrier in the potential energy profile are too weak to yield a concentration profile consistent with experiments. It is interesting to note that the channel is only occupied ~10% of the time, and almost never doubly occupied. The PNP calculations also predict featureless profiles similar to BD, except that cation concentration in the GA channel is enhanced by an order of magnitude compared to the bath solution, which is a direct consequence of the energy well created by neglecting ion self-energies as discussed earlier (Kurnikova et al., 1999; Cardenas et al., 2000).

View larger version (20K): [in this window] [in a new window]   FIGURE 7   Concentration profiles for potassium ions in the GA channel with no applied potential (A) and 200 mV applied potential (B). In both cases, a 500 mM KCl solution is used in the reservoirs. The pore region is divided into 16 equal segments, and each reservoir into two segments. The reservoir concentrations are represented by shaded bars.

In Fig. 8 we show the concentration profile for Ca²⁺ ions obtained under similar circumstances (500 mM CaCl₂ solution, 200 mV applied potential). As emphasized above, the large barrier (Fig. 4 B) prevents the entry of Ca²⁺ ions into the GA channel. Because the barrier gets higher with decreasing _c, we will not see binding sites for any value of dielectric constant. The lack of binding sites for Ca²⁺ ions in GA thus provides the most direct evidence for the failure of continuum electrostatics.

View larger version (15K): [in this window] [in a new window]   FIGURE 8   Concentration profiles for calcium ions in the GA channel, obtained using a 500 mM CaCl2 solution and 200 mV applied potential.

A final piece of evidence demonstrating the failure of continuum electrostatics in the GA channel is its inability to describe the observed saturation of conductance with increasing concentration. This saturation is a direct result of rate-limiting barriers in the channel and will not occur unless the ion has to climb substantial energy barriers with heights >1.5 kT to move out of the energy wells in the channel (Chung et al., 1998). In this respect, the barrier in the GA channel is too small to induce the saturation behavior (Fig. 2 A). A BD study of the current-concentration curve for symmetric KCl solution confirms this expectation, as shown in Fig. 9. The potassium current keeps increasing with concentration with no sign of saturation, in disagreement with the experimental data in which currents clearly saturate with a half-maximum value of K_s  0.23 M (O. S. Andersen, private communication).

View larger version (12K): [in this window] [in a new window]   FIGURE 9   Conductance-concentration curve for potassium ions obtained with a 200 mV applied potential fitted by the solid line.

Although this model accurately reproduces the observed I-V data, it fails to predict the observed binding sites or the saturation of current with increasing concentration. Given these results, one clearly has to look beyond the I-V curves to test the validity of a theoretical approach. On the basis of the complete evidence presented above, we have to conclude that this continuum electrostatics model fails in the GA channel.

Potential energy profiles for GA

The failure of continuum electrostatics in the GA channel leaves MD as the only reliable method for the calculation of forces on ions in this channel. Unfortunately, as noted in the Introduction, none of the MD simulations of GA carried out so far has yielded free-energy profiles that can reproduce experimental currents. The force fields currently used in MD can be improved by including polarization effects, and hopefully, once this is done, MD calculations of potentials of mean force will give more satisfactory results. In the mean time, one can pursue the study of structure-function relationships in GA indirectly by "guessing" the potential profile that will reproduce the available data. Here we give an example of this inverse method by constructing a potential energy profile for potassium ions in the channel.

The fact that GA is a very narrow, single-ion channel (except at high concentrations) makes this task relatively easy. As a first approximation, we can find the profile along the channel axis and supplement it with a harmonic constraint in the radial direction, thus reducing the search to a one-dimensional problem. The shape of the axial profile is, of course, well known from both the experiments and MD simulations. As shown in Fig. 10, it has two binding sites at ~±11 Å, and a central barrier between them. Here we ignore finer details such as small oscillations in the barrier that are not likely to significantly influence the overall conductance properties of GA. Note also that the exact location of the energy wells (for example, placing them at ±9 Å) does not have much bearing on the channel conductance found in BD. This is much more dependent on the depth of the wells, U_w, and the height of the central barrier, U_b. These are two parameters that need to be determined from the BD simulations by fitting the calculated conductance under different applied potentials and concentrations to the available physiological data. In the profile illustrated in Fig. 10, these two parameters, U_w and U_b, are fixed at 8 kT and 5 kT.

View larger version (16K): [in this window] [in a new window]   FIGURE 10   Shape of the potential profile used in the inverse method. The two basic parameters are the depth of the wells Uw and the height of the barrier Ub as indicated by the arrows. Note that the well depth refers to the zero potential in the reservoir but the barrier height is measured with respect to the well minimum. The curved parts are produced using a Fermi function of the form 1/{1 + exp[±(z  z0)/d]}. The width of the well at half-maximum is ~10 Å. The profile shown with Uw = 8 kT and Ub = 5 kT gives the best description of the physiological data on GA.

The diffusion coefficient of K⁺ ions in the GA channel is expected to be considerably suppressed in the GA channel (Roux and Karplus, 1991). However, the spread in the estimated values of D is too large to be able to choose a particular value. Therefore, we consider D as a third parameter to be determined from the BD simulations. The variation of current with D is illustrated in Fig. 11 for the profile with U_w = 8 kT and U_b = 5 kT. The current increases linearly with D at first, but deviates from it with increasing D. This nonlinearity arises from the fact that ions do not diffuse in a flat energy landscape, but have to surmount energy barriers. The results in Fig. 11 demonstrate that the channel conductance can be easily fitted by adjusting the diffusion coefficient of K⁺ ions.

View larger version (12K): [in this window] [in a new window]   FIGURE 11   Variation of the channel currents with the diffusion coefficient of K+ in the channel. The depth of the well and height of the barrier used are, respectively, 8 kT and 5 kT. The results are obtained using an applied potential of 200 mV and an ionic concentration of 500 mM.

To see how the U_w and U_b values influence the channel conductance, we have carried out BD simulations using a wide range of values (Fig. 12). Here an applied potential of 200 mV is used with a 500 mM KCl solution, and D is reduced to 0.1 times the bulk value. Each curve in the figure shows how the channel current decreases exponentially with the increasing barrier height U_b for a fixed value of U_w. Note that when U_b is fixed, the current increases with U_w because a deeper well is more successful in attracting the K⁺ ions into the channel and thereby facilitating their conduction. The broken horizontal line across the figure indicates the experimentally measured current of 4 pA (O. S. Andersen, private communication). Thus, for a given well depth, one can find a matching barrier height that will reproduce the experimental conductance. The requirement that the wells provide binding sites excludes the very low values for U_w (i.e., U_w > 4 kT), but otherwise it does not help in constraining the potential parameters further. This exercise exposes the dangers of relying exclusively on fitting linear I-V curves in model studies. Not only one can fit the channel conductance for a given potential profile by adjusting the diffusion coefficient, but even when D is fixed, there are many sets of potential profiles that can fit the observed conductance. Thus conductance alone cannot provide an adequate test for a phenomenological permeation model.

View larger version (17K): [in this window] [in a new window]   FIGURE 12   Variation of the channel current with the barrier height Ub for fixed values of Uw (indicated at the top of each curve). An applied potential of 200 mV and a 500 mM KCl solution are used in the BD simulations. The diffusion coefficient of potassium ions in the channel is suppressed to 0.1 times the bulk value. Each set of data is fitted by a exponential decay function (solid and dashed lines). The horizontal dashed line indicates the experimental current value from O. S. Andersen (private communication).

We next consider the current-concentration curve and study its sensitivity to the potential parameters and diffusion coefficient. The saturation seen in experimental current-concentration curves cannot be reproduced for all values of diffusion coefficient. For a given D a potential profile cannot necessarily be found that can reproduce the observed saturation. A brief survey of possible diffusion values showed that saturation did not arise when D was larger than 0.3 times the bulk value. This indicates that the diffusion coefficient is more than just a scaling factor, it plays a dynamical role in ion permeation, influencing the saturation properties of the channel.

A successful description of saturation, however, is possible when the potential parameters are chosen as U_w = 8 kT, U_b = 5 kT, and D is reduced to 0.1 times the bulk value. In the remaining figures, we study in greater detail the consequences of this chosen parameter set. In Fig. 13, the saturation curves obtained under two different applied potentials (100 and 200 mV) are compared to the experimental data of O. S. Andersen (private communication). The calculated curves follow the experimental result closely at both driving potentials. The simulation results are fitted with a Michaelis-Menten curve as indicated by the solid lines in the figure. The calculated half-saturation value of 250 mM is in good agreement with the experimental value of 230 mM at 200 mV. The corresponding values at 100 mV are 90 mM (calculated) and 105 mM (experimental). Fig. 14 shows an I-V curve obtained from simulation of a symmetric 500 mM KCl solution. The results of our simulations (filled circles) are compared with the experimental measurements (open circles) obtained by O. S. Andersen (private communication). The agreement with the data is again very good.

View larger version (14K): [in this window] [in a new window]   FIGURE 13   Saturation of current with concentration at two applied potentials, 100 mV (bottom curve) and 200 mV (top). The profile in Fig. 10 is used in BD simulations with Uw = 8 kT, Ub = 5 kT, and D reduced to 10% of the bulk value. The experimental data of O. S. Andersen (private communication) are superimposed on the figures (open symbols), and are fitted by dashed lines..

View larger version (11K): [in this window] [in a new window]   FIGURE 14   Current-voltage relationship obtained with BD simulations using the profile in Fig. 10. The results obtained by using 500 mM KCl solution (black circles fitted with a line) are compared with the experimental data of O. S. Andersen (private communication), shown by open circles.

Finally, in Fig. 15, we show the concentration profiles for K⁺ ions in the GA channel obtained from a symmetric 500 mM KCl solution. As expected, the ion concentration is very large at the binding sites and depressed in the middle, regardless of whether there is an applied potential of 200 mV. The magnitude of the concentration may appear too large at first sight. This is simply due to the small volume available at the binding sites. In fact, there are only ~0.75 ions on average at each site in the case shown. This also indicates that the channel is often multiply occupied, unlike when the flatter profile found from continuum electrostatics is used, as the larger wells more easily attract ions into the channel and retain them at the binding sites. The number of ions in the channel increases with concentration, as at higher concentrations ions can find their way into the channel more quickly. For example, at 150 mM there are on average 0.8 ions in the channel, while at 1 M there are 1.75. An analysis of ion trajectories in the multiply occupied channel shows that the inter-ion separation remains roughly constant for any ion pair.

View larger version (22K): [in this window] [in a new window]   FIGURE 15   Concentration profiles for potassium ions in the GA channel as in Fig. 7, but using the profile in Fig. 10, with no applied potential.

While it is impossible to prove the uniqueness of the parameter set we have obtained, our results nevertheless suggest that large variations from these values are unlikely to lead to a satisfactory description of the data. For example, Fig. 12 indicates that the potential parameters U_w = 12 kT, U_b = 5 kT may also work, although the larger wells will increase the probability of multiple occupancy. It is also of interest to compare our chosen profile with those of Chiu and Jakobsson (1989). They reproduce sodium conductance properties via electrodiffusion equations using a similar potential profile with the parameters U_w = 5.4 kT, U_b = 4.2 kT, and a diffusion coefficient for Na⁺ ions that is 0.07 times the bulk value. McGill and Schumaker (1996) also find that similar well depths and barrier heights are required to match experimental conductances using their diffusion theory. Thus, there is a reasonable congruence among all the sets of parameters.

The saturation property of the GA channel is seen to provide the most sensitive quantity for the purpose of determining the potential energy profile of ions. In most model studies of ion channels, conductance and I-V curves are studied in great detail while no attention is paid to the saturation curve. In fact, as illustrated here, reproducing the saturation curve provides a more stringent test for a permeation model and should be given more consideration in future studies.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Our main result is that continuum electrostatics using a rigid protein structure cannot provide a consistent description of all the available data on the GA channel. Surprisingly, the use of a dielectric constant of 80 in the channel appears to give the best results in our continuum electrostatics, despite the fact that from microscopic considerations such a high value seems unreasonable in the narrow GA channel. Although slight changes in the protein structure may alter the energy profiles, the experimental data limit the amount of flexibility allowed, and it seems unlikely that all the problems encountered could be overcome using a different or flexible protein. Given the obvious problems in describing polarization around an ion and the effect of water dipoles, we believe that continuum electrostatics should not be used to model GA. This applies to any model that relies on the solution of Poisson's equation using a dielectric continuum, such as the current BD simulations (in which the forces are calculated from Poisson's equation) and continuum theories such as Poisson-Boltzmann and Poisson-Nernst-Planck equations. This leaves MD as the only method for calculating the forces on ions in the GA channel. Unfortunately, the force fields currently used in MD studies are too crude for this purpose and require further refinement to be able to make contact with experimental data. We have used the inverse method to construct a potential energy profile for potassium ions that gives an adequate description of the available physiological data. It will be interesting to compare this profile with future profiles to be determined from MD with improved force fields. In the meantime, the inverse approach can be utilized to relate the structure-function data in the GA channel. Studies of asymmetric and bi-ionic solutions could provide further tests of our proposed potassium potential profile. Also, we expect that profiles for other monovalent ion species can be constructed to match their conductance values by making small changes to the potassium profile. Work in this direction is in progress.

An important question is whether the above results for the GA channel have any bearing on biological ion channels. Because GA was the only channel with a well-defined tertiary structure for a long time, it became a model case for all ion channels. The hope was that insights gathered from the study of the GA channel could be used in understanding the permeation properties of other channels. The determination of the KcsA structure (Doyle et al., 1998) and subsequent studies of the permeation mechanisms in potassium (Chung et al., 1999, 2002) and calcium channels (Corry et al., 2001) reveal that GA is not likely to fulfill this role. Both in terms of structure (single filing of water versus presence of water-filled cavities and vestibules) and ion dynamics (single-ion versus multi-ion permeation mechanism), GA is very different from these biological channels. It is the long single-file chain of water molecules that creates problems when using continuum electrostatics to model the GA channel, as the use of uniform dielectric constants cannot mimic their long-range polarization. In contrast, the continuum electrostatics-BD approach has been used successfully to reproduce a wide variety of experimental data when modeling biological ion channels, including the KcsA potassium channel (Chung et al., 1999, 2002). The reason for this is twofold: first, continuum electrostatics works reasonably well in the wider cavity and vestibule regions that form the major part of these channels; second, the selectivity filter regions where the single filing occurs are much shorter, involving only a few water molecules sandwiched between ions. That is, there are no long chains of water molecules as in the GA pore that appear to cause the problems seen here. Thus, we believe that GA is a very special channel that needs to be handled with extreme care. The rewards from its study are not to be found in getting direct insights about biological ion channels, but rather constructing reliable models of permeation in a difficult test case that can later be applied to other channels.