INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

A central problem in studies of ion permeation through biological membrane channels is to understand how channels can be both highly selective and yet still conduct millions of ions per second. Calcium channels exemplify this problem; they are ubiquitous in excitable cells and extremely selective, selecting calcium over sodium at a ratio of 1000:1 (Hess et al., 1986); yet the picoampere currents they support require >10⁶ calcium ions to pass per second (Tsien et al., 1987). Unlike potassium channels, which have a narrow selectivity filter and only allow ions of a particular size to pass (Doyle et al., 1998; Allen et al., 1999a, 2000a), calcium channels select between ions of almost identical radius, the Pauling radii of sodium and calcium ions being 0.95 and 0.99 Å, respectively. Moreover, calcium channels are known to admit much larger ions, the largest observed is tetramethylammonium, with a radius of ~2.8 Å (McCleskey and Almers, 1985). Thus, a different mechanism of selectivity from that in the potassium channel must be at play, one that relies on the different charges on the ions. Monovalent ions can permeate the channel in the absence of calcium at much higher levels of conductance than can any divalent ions (Kostyuk et al., 1983; Almers and McCleskey, 1984; Hess et al., 1986; Kuo and Hess, 1993a), but are blocked when the calcium concentration reaches only 1 µM (Kostyuk et al., 1983; Almers et al., 1984). That this block is dependent on membrane voltage (Fukushima and Hagiwara, 1985; Lansman et al., 1986; Lux et al., 1990) and the direction of ion movement (Kuo and Hess, 1993a, b) has been taken as evidence for a multi-ion binding (or selectivity and blocking) site residing in the pore. Four glutamate residues in close proximity are believed to line the pore and to be a component of the selectivity filter of the channel, as point mutations of these change the characteristics of selectivity (Yang et al., 1993; Kim et al., 1993; Ellinor et al., 1995; Parent and Gopalakrishnan, 1995; Bahinski et al., 1997). The glutamate residues are expected to be highly charged and to strongly bind the calcium ions in the channel leading them to block the passage of sodium ions.

A number of theoretical models have been developed to explain permeation and selectivity in the calcium channel. Single-file rate theory models in which ions sequentially hop from one site to another have been used most extensively (Tsien et al., 1987). Because of the difficulty in obtaining both high selectivity and throughput with a single binding site (Bezanilla and Armstrong, 1972), these models originally contained two sites in which repulsion between ions in neighboring sites increases transit rates (Hess and Tsien, 1984; Almers and McCleskey, 1984). As the two-site models could not accommodate the mutation data, a new rate model was recently proposed where a single-site is flanked by lower affinity sites to aid the exit of ions from the central site (Dang and McCleskey, 1998). Other mechanisms involving single sites have also been developed, such as competition between calcium ions for the binding charges (Armstrong and Neyton, 1991; Yang et al., 1993). These rate theory models have provided many useful insights as to how calcium channels may achieve their selectivity with a high throughput. However, they cannot be used to relate the structural parameters of the channel to functional elements (McCleskey, 1999). For example, in these theories no physical distances or shapes are used and there is no direct connection between energy minima used in the theory and physical sites in the pore.

A first attempt to relate the observed properties of the calcium channel to its structure was made with the Poisson-Nernst-Planck (PNP) theory, which uses continuum electrostatics and electrodiffusion equations to calculate channel conductance (Nonner and Eisenberg, 1998). The shortcomings of the PNP theory as applied to a model calcium channel were pointed out by McCleskey (1999) and Miller (1999). These criticisms have been given a solid foundation in recent comparisons of PNP theory with Brownian dynamics (BD) simulations (Corry et al., 1999, 2000a, b), which show that the mean field approximation used in the PNP theory completely breaks down in narrow channels, such as the calcium channel. The good agreement between the PNP results and the channel data, often put forward as a proof of its validity, is seen in hindsight as a fortuitous outcome of mixing incorrect physics with unrealistic parameter values. For example, the calcium diffusion coefficient used in the PNP fits (10⁵ times the bulk value) is 10,000 times smaller than the microscopic estimates obtained from molecular dynamics simulations, which suggest at most a 10-fold reduction in calcium diffusion compared to the bulk value (Allen et al., 2000b). Agreement with experiment also relies on the inclusion of ad hoc chemical potentials whose electrostatic origin is not clear.

Failure of the mean field approximation in narrow channels indicates that any theory that aspires to relate channel structure to its function must treat ions explicitly. Because all the atoms in the system are treated explicitly in molecular dynamics, it would provide the ultimate approach to the structure-function problem. Unfortunately, computation of most channel properties (e.g., conductance) using molecular dynamics is still beyond the capabilities of current computers. The only remaining alternative is BD simulations, where water is treated as continuum and only the motion of individual ions is followed via the Langevin equation. The early BD simulations of ion permeation were carried out in one dimension with assumed electric potentials (Cooper et al., 1985; Bek and Jakobsson, 1994), which were not very useful as realistic models of channels. In the past few years we have extended BD simulations to three dimensions, with the electric fields properly calculated from the solution of Poisson's equation (Li et al., 1998; Hoyles et al., 1998a). These realistic BD simulations have been used to describe ion permeation in the acetylcholine receptor (Chung et al., 1998) and KcsA potassium channel (Chung et al., 1999). Multi-ion interactions were found to be instrumental in explaining the high throughput of potassium channels, and are expected to play a similarly significant role in understanding the high conductance of calcium channels.

The aim of this paper is to construct a simple model of the structure of calcium channels and examine its various properties using electrostatic calculations and BD simulations. The parameters in the model are determined from either molecular dynamics or a variational principle that optimizes the quantity in question. Thus there are no free parameters that are fitted to data, nor ad hoc chemical potentials that are arbitrarily chosen. The model relates structural features to functional roles and, as will be seen, successfully predicts many of the observed properties of the calcium channel using only the principles of electrodynamics.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Channel model

The crystal structure of calcium channels is not known at present. Nevertheless, through a judicious use of important clues from various experiments one can develop a simplified model of the calcium channel that should be sufficiently accurate for the purposes of electrostatic calculations and BD simulations. The cross-section of the channel model used in this work is shown in Fig. 1 A. A three-dimensional shape of the channel is generated by rotating the curves in Fig. 1 A about the axis of symmetry (z axis) by 180°. The channel extends from z = 25 Å to 25 Å, long enough to span a typical membrane. In constructing this model, we have followed the basic topology of the potassium channel (Doyle et al., 1998); that is, a narrow selectivity filter, connected to a wide chamber that tapers off as it approaches to the intracellular side. One significant difference from the potassium channel is the existence of a relatively short vestibule on the extracellular side with a fairly wide opening. This is suggested by molecular modeling studies (Schetz and Anderson, 1993; Guy and Durell, 1995; Doughty et al., 1995, 1998) of the known amino acid sequences of the calcium channel (Tanabe et al., 1987; Mikami et al., 1989; Williams et al., 1992). A larger external mouth compared to the internal one is required to explain the asymmetry between the inner and outer saturation currents (Kuo and Hess, 1992).

View larger version (18K): [in this window] [in a new window]   FIGURE 1   Model calcium channel and reservoirs. (A) The three-dimensional channel model is generated by rotating the curves about the central axis by 180°. The positions of two of the four glutamate groups are shown by the squares, and the inner end of two of the four mouth dipoles by the diamonds. The other two groups lie into and out of the page. The intracellular end of the channel is on the left and the extracellular side on the right. (B) The channel is enclosed with cylindrical reservoirs on either side representing the intracellular and extracellular baths.

The radius of the selectivity filter is determined from the size of the largest permeable ion (tetramethylammonium) as 2.8 Å (McCleskey and Almers, 1985). Interpretation of the mutation data in reaction rate theories suggests that the four glutamate residues (EEEE locus) in the selectivity filter must be in close proximity in order to form a single binding site (Yang et al., 1993; Ellinor et al., 1995; Bahinski et al., 1997). This is further supported by the voltage dependence of calcium block, which suggests that calcium binds at the same location whether entering the channel from the inside or outside (Kuo and Hess, 1993a). Therefore, we have chosen the length of the selectivity filter to be 5 Å, which is much shorter compared to that in a potassium channel (12 Å). The position of the selectivity filter in the channel is not known, though it is suspected to be toward the external side of the channel, as it is more accessible to ions from the outside of the channel than from the inside (Kuo and Hess, 1993a). Our trials with various positions of the selectivity filter in the channel also confirm this conjecture: when the filter position is further removed from the external mouth, it is not possible to reproduce most of the known properties of calcium channels. The wide chamber near the middle plays a similar stabilizing role to that in potassium channels, providing a water-filled cavity for ions exiting from the selectivity filter (Roux and Mac-Kinnon, 1999).

The highly charged glutamate residues forming the selectivity filter play an essential role in determining the channel conductivity and selectivity, and therefore, choosing their positions and charges correctly is of critical importance. The four glutamate residues are modeled by four fixed charges located in close proximity, but spread asymmetrically in a spiral pattern 1 Å behind the channel wall. The placement of charges in an asymmetric pattern rather than in a ring helps to account for the mutagenesis studies that show the removal of each charge has a different effect on channel conductance. The four charges are located at z = 10.50, 11.83, 13.17, and 14.5 Å, and each rotated by 90° from the last (only two are shown in Fig. 1 A). Finally, to overcome the large image forces at the intracellular end of the channel we have placed four mouth dipoles, 5 Å in length, with their inner ends 1 Å inside the pore wall at z = 17.5 Å. The charges on glutamates and mouth dipoles are optimized to obtain the maximum ionic currents as discussed below. Because we use a rigid protein structure, we do not consider here the possibility of negative charges protruding into the channel and swinging out of the way as ions pass through. We find that such a flexibility of glutamate residues is not required to reproduce experimental data.

The dielectric constant of the channel protein is taken uniformly as  = 2. The dielectric constant of water inside the channel environment is not well known, as it is difficult to determine its value directly from experiments. Recent molecular dynamics simulations of water inside narrow channels have suggested that it may be considerably lower than its bulk value (Sansom et al., 1997). However, BD simulations of ion permeation in potassium channels indicate that current ceases to flow if  in the channel is lower than 40 (Chung et al., 1999). In view of these uncertainties, we have adopted the value of  = 60 that allows large conductance through the model channel. Further justification for this choice will be given later.

Reservoirs

The channel in Fig. 1 A is enclosed by cylindrical reservoirs on either side of the membrane that represent the intracellular and extracellular baths (Fig. 1 B). The reservoirs have a fixed radius of 30 Å and their height is adjusted to obtain the desired concentration (typically ~33 Å). This length scale (4 Debye lengths for 150 mM) is optimal in the sense that the fields of ions outside a 30 Å sphere are totally screened out, and therefore they would have no effect on the dynamics of ions inside the channel. The reservoir boundaries simply serve to confine the ions within the simulation system, which is the easiest way to maintain the average concentrations in the baths at the desired values. Implicit in the use of reservoirs is the assumption that electrolyte solution continues beyond its boundaries. It is worthwhile to emphasize that the reservoir boundaries are not used in the solution of Poisson's equation. That is, we do not fix the potentials at the top and bottom reservoir surfaces, rather they follow from the solution of Poisson's equation as described below.

Scattering of an ion from the boundary wall can be viewed as an ion moving out of the reservoir and another one entering at the same time. In reality, the number of ions in the reservoir will not be constant, but fluctuate around an average value. Such fluctuations have recently been taken into account using a grand canonical Monte Carlo method (Im et al., 2000). Here we use the simpler method because fluctuations in concentration are expected to increase the noise in the current measured from BD simulations but not to affect its average value, which is the observed quantity.

Solution of Poisson's equation

The electric potential  of a configuration of ions and fixed charges in the channel system, represented by the charge density , is found from the solution of Poisson's equation:

(1)

where the dielectric constant (r) has different values on either side of the channel boundary. For the proposed boundary in Fig. 1 A, Poisson's equation can only be solved numerically. This is achieved using the boundary charge method (Levitt, 1978), where the boundary is divided into small sectors, and each sector is represented by a point charge at its center. For faster convergence and more accurate solutions, we have included the effect of curvature of sectors in the solutions following Hoyles et al. (1996, 1998b). We refer to these references for details of this method.

Implementing different values of  for water in the channel and reservoirs leads to problems in the boundary charge method (Chung et al., 1999). We use instead the same value of  = 60 for water everywhere and incorporate the neglected Born energy difference between the reservoirs and the channel interior as a potential barrier at either channel entrance. The barrier height is estimated from the solution of Poisson's equation using a three-dimensional grid as described elsewhere (Moy et al., 2000). Further description and justification of this approximate way of handling the dielectric constant in BD simulations is given in Chung et al. (1999). We note that while solution of Poisson's equation using a three-dimensional grid avoids the above problem, it is far too slow to be of practical use in BD simulations.

Brownian dynamics

An introduction to BD simulations in one-dimensional channels is given by Cooper et al. (1985). BD simulations have recently been extended to realistic three-dimensional channel geometries (Li et al., 1998; Chung et al., 1998, 1999; Hoyles et al., 1998a). We give a brief description of the method here and focus on new features that have not been discussed before.

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

(2)

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. 2 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 (Hoyles et al., 1998a). For this purpose, the total electric potential _i experienced by an ion i is broken into four pieces

(3)

where the sum over j runs over all the other ions in the system; _X,i is the external potential due to the applied field, fixed charges in the protein wall, and charges induced by these; _S,i is the self-potential due to the surface charges induced by the ion i on the channel boundary; _I,ij is the image potential due to the charges induced by the ion j; and _C,ij is the Coulomb potential due to the ion j. The first three potential terms in Eq. 3 are stored in, respectively, 3-, 2-, and 5-dimensional tables (dimension is reduced by one in the latter two cases by exploiting the azimuthal symmetry of the dielectric boundary).

The short-range forces are used to keep the ions in the system and also to mimic other interactions between two ions that are not included in the simple Coulomb interaction. In order to prevent ions from leaving the system, a hard-wall potential is activated when the ions are within one ionic radius of the reservoir boundaries, which elastically scatters them. For the ion-wall interaction U_IW, we use the usual 1/r⁹ repulsive potential

(4)

where R_i is the ion's radius, R_w is the radius of the atoms making up the wall, R_c(z) is the channel's radius as a function of the z coordinate, and a is the ion's distance from the z axis. We use R_w = 1.4 Å and F₀ = 2 × 10¹⁰ N in Eq. 4, which is estimated from the ST2 water model used in molecular dynamics (Stillinger and Rahman, 1974).

At short ranges, the Coulomb interaction between two ions is modified by adding a potential U_SR(r), which replicates effects of the overlap of electron clouds and hydration. Molecular dynamics simulations show that the hydration forces between two ions add further structure to the 1/r⁹ repulsive potential due to the overlap of electron clouds in the form of damped oscillations (Guàrdia et al., 1991a, b). These two effects can be approximately represented by

(5)

Here the oscillation length c_w = 2.76 Å is given by the water diameter and the other parameters are determined by fitting Eq. 5 to the potentials of mean force given by Guàrdia et al. (1991a, b, 1993, 1996). Fig. 2 A shows a plot of the short-range potential for NaCl solution used in our Brownian dynamics simulations. For anion-cation pairs R_c = r₁ + r₂, but for like ions the contact distance is pushed further to R_c = r₁ + r₂ + 1.6 Å. The origin of the hydration force R is slightly shifted from R_c; by +0.2 Å for like ions and by 0.2 Å otherwise. The exponential drop parameter is determined as c_e = 1 Å for all ion pairs. Finally, the overall strength of the potential is U₀ = 16.8, 8.5, 1.7, 2.5, 0.8, and 1.4 kT for Ca-Cl, Na-Cl, Ca-Na, Na-Na, Ca-Ca, and Cl-Cl pairs, respectively. This potential agrees well with the potential of mean force derived by Guàrdia et al. (1991a, b). The short-range force in Eq. 2 is determined from the derivatives of the potentials in Eqs. 4 and 5.

View larger version (18K): [in this window] [in a new window]   FIGURE 2   Ion-ion forces used in BD simulations. (A) The inter-ion potentials for Na+-Na+ (solid line), Na+-Cl (dashed line), and Cl-Cl (dotted line) ion pairs are plotted against the ion separation as given by Eq. 5. (B) The radial distribution functions for 1.79 M NaCl solution derived from BD simulations (the same line styles as in A are used). The locations of the maxima found in the molecular dynamics simulations of Lyubartsev and Laaksonen (1995) are indicated by the arrows at the top of the graph.

The BD simulations using this combination of Coulomb and short-range forces accurately mimic the results of molecular dynamics simulations. In Fig. 2 B we show the radial distribution functions for Na-Na, Na-Cl, and Cl-Cl pairs obtained from a 2.5 ns BD simulation of a 1.79 M NaCl solution (22 Na⁺ and 22 Cl ions) confined in a large cylinder with a diameter and height of 30 Å. To avoid the edge effects, ions within 8 Å of the boundary at any time step are excluded from the sampling. As expected, the resulting peaks in the distribution function are located at the minima of the potential of mean force, and also closely match those locations found in the radial distribution functions from molecular dynamics simulations of Lyubartsev and Laaksonen (1995) (indicated by the arrows in Fig. 2 B). Similar results are obtained for CaCl₂ solutions, but because the molecular dynamics simulation results obtained by different groups vary, a comparison similar to Fig. 2 B is not presented for CaCl₂.

We have found that simpler ion-ion interactions used in BD studies of other channels (e.g., Chung et al., 1998, 1999) are not suitable for use in calcium channels. Contrary to the realistic interaction described above, they allow cations to pass each other in the selectivity filter, thus making it impossible to explain the observed blocking of sodium ions by calcium, and vice versa.

The Langevin equation (Eq. 2) is solved at discrete time steps following the 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 (between z = 25 to 20 Å) 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 also assigned randomly 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. For this purpose, the ion on the furthermost right-hand side is chosen, and it is placed on the far left-hand side of the left reservoir, making sure that it does not overlap with another ion. The stochastic boundary trigger points, located at either pore entrance, are checked at each time step of the simulation. The sudden disappearance of an ion from the reservoir boundary has a negligible effect on ions in the channel. To give an example, the force between two monovalent ions at 30 Å is 3 × 10¹³ N, which is 1000 times smaller than the average random force. The potential change experienced by an ion near the trigger point in the channel is about 3 mV when an ion is removed from one reservoir and placed in the other.

The BD program is written in FORTRAN, vectorized and executed on a supercomputer (Fujitsu VPP-300). 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 roughly 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, , in Eq. 2 by the Einstein relation.)

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Channel parameters

The three channel parameters that are not known experimentally and need to be determined by other methods are the magnitude of the charges on glutamate residues and mouth dipoles, the dielectric constant of water, and the diffusion coefficient of ions in the channel. A straightforward fit of these parameters to the available data is not very satisfactory, since one is likely to find many possible sets that eventually have to be distinguished on their physical merits. Therefore, we prefer using guiding principles such as optimization or a more explicit theory (e.g., molecular dynamics) in estimating these quantities.

The determination of the molecular structure of the proteins may help to find the magnitude of the charge of residues in the channel. In the meantime, we expect that the charges in the channel would have evolved to maximize the transit rate of calcium ions. In Fig. 3 A we show the dependence of the calcium current on glutamate charges. The BD simulation results in this figure are obtained using symmetric 150 mM CaCl₂ or NaCl solutions with an applied field of 2 × 10⁷ V/m (corresponding to a potential of ~200 mV producing an inward current). As the charge on the glutamate groups is systematically increased (while the charge on the mouth dipoles is held fixed), the calcium current found from BD simulations sharply increases from zero to a narrow peak at a charge of 1.3 × 10¹⁹ C before dropping steeply again at greater charge strengths. In fact, no calcium current is measured during our simulations if the charges are <1.0 × 10¹⁹ C or >1.6 × 10¹⁹ C. The sodium current also peaks at the same value but conducts over a greater range of glutamate charges, as is shown by the open circles. It is noteworthy that the peak calcium current occurs for such a narrow range of glutamate charges. A fully charged glutamate group has a charge of e (1.6 × 10¹⁹ C). However, in an electrolyte solution the charges are likely to become protonated, leading to a lower effective charge on the residues. (Chen et al., 1996; Morrill and MacKinnon, 1999; Root and MacKinnon, 1994, Chen and Tsien, 1997). As the amount of protonation is not known, we use the optimum value of the glutamate strengths, 1.3 × 10¹⁹ C, for the remainder of this study.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Dependence of channel current on fixed charge strengths. (A) The current passing through the channel with 150 mM CaCl2 (filled circles, left side scale) and 150 mM NaCl (open circles, right side scale) under a 200 mV driving potential is plotted against the charge on each of the glutamate groups. The magnitude of the charge on the mouth dipoles is fixed at 0.6 ×1019 C. Filled circles are obtained from a 1.0 µs and open circles from a 0.5 µs simulation period. (B) The outward (filled circles) and inward (open circles) current passing through the channel with 150 mM CaCl2 in the reservoirs and a 200 mV driving force is plotted against the magnitude of the charge on each end of the mouth dipoles. The charge on the glutamates is held at 1.3 ×1019 C. Results are obtained from a 2 µs simulation period. Error bars in this and following figures have a length of one standard error of the mean and are not shown when smaller than the data points.

In a similar investigation, the strength of the mouth dipoles is varied while the glutamate charges are held fixed at their optimal value. As shown in Fig. 3 B, the outward calcium current is critically dependent on the charge strengths, as in the case of the glutamate residues. The current is maximum when a charge of 0.5 × 10¹⁹ C is placed on each of the four dipoles. A further increase in the dipole strength reduces the current rapidly (filled circle, Fig. 3 B). In contrast, the inward current exhibits a different dependence on the mouth dipole strength (open circles, Fig. 3 B). The current increases steeply with the dipole strength and then remains constant with a further increase. In all subsequent simulations, we use a charge of 0.6 × 10¹⁹ C, which falls between the optimum values of inward and outward currents and gives close to the maximum value for each.

Molecular dynamics studies of water in spherical cavities (Zhang et al., 1995) and narrow pores (Sansom et al., 1997) suggest that the dielectric constant  is substantially reduced from the bulk value. The effect of changing the dielectric constant on the results of BD simulations in narrow pores (the potassium channel) has been examined elsewhere (Chung et al., 1999). This study also found that the optimum charge strengths are insensitive to the value of the dielectric constant. The dielectric constant of water in the channel is chosen as  = 60. While this value is rather close to the bulk value, the channel ceases to conduct calcium ions if lower values of  are used. For example, when the dielectric constant  inside the channel is assumed to be 50 and a potential difference of 200 mV is applied, the current across the channel is only 2.4 ± 0.6 pA, compared to 7.1 ± 0.6 pA with  = 60. With a further reduction of  to 40, the current is reduced to 0.4 ± 0.2 pA (during a simulation period of 3 µs). Virtually no conduction takes place with an applied potential of 100 mV and  of 50. In a simulation period of 5.5 µs, only one calcium ion crosses the channel, resulting in a current of 0.06 pA.

The diffusion coefficient of ions inside the channel can be estimated from molecular dynamics simulations. There are a number of such studies which indicate that the diffusion coefficient is significantly reduced from its bulk value inside narrow channels (Roux and Karplus, 1991; Lynden-Bell and Rasaiah, 1996; Smith and Sansom, 1997, 1999; Allen et al., 1999a, b, 2000b). Allen et al. (2000b) have carried out a systematic study of diffusion coefficients of K⁺, Na⁺, Ca²⁺, and Cl ions in cylindrical channels with radii varying from 3 to 7 Å. Here we use their estimates as a guide and use 0.5 times the bulk diffusion coefficient for calcium ions in the channel chamber (25 < z < 7.5 Å) and 0.1 times the bulk in the selectivity filter (7.5 < z < 20 Å). Corresponding values of 0.5 and 0.4 times the bulk value are used for sodium. The bulk values are used in the reservoirs for all ions.

In Fig. 4 we illustrate the sensitivity of the channel conductance on the choice of diffusion coefficient. Here the diffusion coefficient of calcium ions is systematically varied from 0.05 to 0.5 times its bulk value in the selectivity filter while it is kept at 0.5 times the bulk in the chamber, and the resulting current is plotted. Contrary to intuitive expectations from continuum theories, the current does not increase linearly with the diffusion coefficient but rather saturates as one approaches toward the bulk value. For example, at the chosen value of 0.1 times the bulk, the calcium current is suppressed by only a factor of 2 rather than 10. This happens because ions have to overcome a potential barrier to conduct, which causes saturation of current in the limit of high diffusion coefficients (ballistic limit). Thus, we expect the results presented in this paper to be quite robust against variations in the diffusion coefficients.

View larger version (13K): [in this window] [in a new window]   FIGURE 4   Dependence of calcium current on the ion diffusion coefficient in the narrow neck region of the channel (7.5 Å < z < 20 Å) plotted as a fraction of its bulk value (0.79 × 109 m2 s1). A concentration of 150 mM CaCl2 is maintained in the reservoirs and a 200 mV driving force is used. Results are obtained from a 2 µs simulation period.

Permeation of calcium and sodium ions

The ion-channel and ion-ion interactions hold the clue to understanding ion permeation mechanisms in channels. Therefore, we first present a detailed study of multi-ion potential energy profiles in the model channel to gain some useful insights. A quantitative description of ion permeation that includes the effects of the thermal motion of ions and their interaction with water molecules requires a dynamic approach, which will be discussed in the following sections by performing BD simulations.

Energy profiles

View larger version (16K): [in this window] [in a new window]   FIGURE 5   Electrostatic energy profile of an ion traversing the channel. The potential energy of an ion held at 1 Å intervals in the z direction but allowed to move to its minimum energy position in the x and y directions is plotted for a calcium ion (solid line, a) and a sodium ion (dashed line, b) in the absence of any fixed charges. When the glutamate groups and mouth dipoles are included, as shown in the inset, the profiles are replotted for calcium (solid curve, d) and sodium ions (dashed line, c). No applied potential is used. We note that 1 kT = 4.11 × 1021 J.

View larger version (14K): [in this window] [in a new window]   FIGURE 6   Electrostatic energy profiles with two calcium ions in the channel under a 100 mV driving force. The potential energy of a calcium ion entering the channel is calculated at 1 Å intervals along the z axis while another calcium ion is resident in the filter (dashed curve). Similarly, the potential energy encountered by the left-hand calcium ion as it attempts to cross-the channel is calculated at 1 Å intervals (solid curve). The second ion is allowed to move to its minimum energy position in the narrow channel neck in both cases. The equilibrium positions of the two calcium ions in the channel are indicated by the arrows. It should be noted that these are two distinct curves and the driving potential cannot be calculated from the total energy drop from right to left.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   The energy profiles as in Fig. 6 except for a sodium ion with one (lower curves) or two (upper curves) other sodium ions in the channel neck. An applied potential of 100 mV is used. The equilibrium positions of the resident sodium ion(s) when the test ion is in the reservoir are indicated by the upward arrows.

Current-voltage relationships

View larger version (10K): [in this window] [in a new window]   FIGURE 8   Current-voltage relationships. The magnitude of the current passing through the channel with a symmetric solution of (A) 150 mM CaCl2 and (B) 150 mM NaCl in both reservoirs is plotted against the strength of the driving potential. The experimental results of Rosenberg and Chen (1991) in similar conditions are shown in the insets for comparison. A simulation period of 4 to 8 µs is used for calcium and 0.5 µs for sodium.

Ions in the channel

View larger version (17K): [in this window] [in a new window]   FIGURE 9   Average number of ions in the channel with an applied potential of 200 mV. The channel is divided into 30 sections, as shown in the inset, and the average number of ions in each calculated over a simulation period (0.5 µs) with (A) 150 mM CaCl2 and (B) 150 mM NaCl in the reservoirs.

View larger version (17K): [in this window] [in a new window]   FIGURE 10   Average number of ions in the channel as in Fig. 9 except with an applied potential of +200 mV with (A) 150 mM CaCl2 and (B) 150 mM NaCl in the reservoirs.

View larger version (10K): [in this window] [in a new window]   FIGURE 11   Rate-limiting steps for ion permeation. The energy profile presented to a calcium ion as in Fig. 6 and the main time-consuming steps for ion permeation are shown for (A) a 200 mV and (B) a +200 mV driving potential. In A the ions permeate from right to left and meet a central energy barrier VB = 2.9 kT. In B the ions permeate in the opposite direction and meet a barrier of 3.7 kT. The time taken for a second ion to enter the channel, 1, and the time for an ion to cross the central barrier, 2, are indicated.

Conductance-concentration relationships

(6)

View larger version (14K): [in this window] [in a new window]   FIGURE 12   Current-concentration relationships. The current obtained with symmetrical solutions of varying concentrations of (A) CaCl2 (filled circles) and (B) NaCl (open circles) in the reservoirs. An applied field of 200 mV is used and the data points are fitted by the solid line using Eq. 9. In A the experimental data of Hess et al. (1986) are shown by the open diamonds and dashed line for comparison. Note that the different scales on the simulation and experimental results are largely due to the different applied potentials in each case. For the BD results a simulation period of 4-8 µs and 0.5 µs are used for calcium and sodium, respectively.

Mixtures of calcium and sodium ions

It is important to see whether our model channel can account for experimental results with more than one ion species present. In particular we look at mixtures of calcium and sodium ions as an example of selectivity between monovalent and divalent ions. To answer such questions as the effect of each type of ion on the permeation of the other and the competition between different types of ions to access the selectivity filter, we again first consider potential energy profiles for mixed ions and then carry out BD simulations.

Energy profiles for mixed ions

View larger version (13K): [in this window] [in a new window]   FIGURE 13   Energy profiles indicating calcium block. The right curve (dashed line) shows the potential energy of a sodium ion given that there is a calcium ion in the filter as indicated in the inset. The left curve (solid line) shows the potential energy of a calcium ion given that there is a sodium ion in the filter. The energies are calculated at 1 Å intervals as in Fig. 5 under a 100 mV driving potential.

View larger version (12K): [in this window] [in a new window]   FIGURE 14   Energy profiles as in Fig. 13 except with a calcium ion on the right side of either (A) one, or (B) two sodium ions as shown in the inset.

Current-voltage relationship in mixed solutions