INTRODUCTION

Top Abstract Introduction Methods Results Discussion References

Models of ion transport in membrane channels have been dominated by two competing approaches based on macroscopic continuum and reaction-rate theories (Hille, 1992). The first involves the solution of two coupled differential equations, the Nernst-Planck equation for the diffusion of ions and Poisson's equation for the mean electric potential. In the second approach, popularly known as Eyring rate theory, the diffusion of ions is modeled as a hopping process. That is, ions move by jumping from site to site, the probability of a hop being determined by the energy difference between the sites and the available thermal energy. Both approaches have their shortcomings when applied to ion channels, and we refer to Cooper et al. (1985) for a critical review. Cooper et al. (1985) made a case for use of Brownian dynamics (BD) in describing ionic motion in membrane channels. In BD, one simulates the motion of ions by solving the Langevin equation. Collisions of an ion with the surrounding water molecules are mimicked by a random fluctuating force plus an average frictional force. The systematic electric forces acting on ions due to various sources can be determined, in principle, by solving Poisson's equation for the system of ions and the protein boundary. Then by iterating the Langevin equation at small time steps, it is possible to follow the trajectories of ions. Thus, as long as one can neglect the dynamics of water molecules in a channel, BD provides a direct solution for the many-body problem of ions confined within an arbitrary boundary.

Despite the strong plea made by Cooper et al. (1985), there had been only a few BD studies of ion channels, mainly by Jakobsson and collaborators (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and Jakobsson, 1994). In all of these studies, model channels were taken as narrow cylindrical tubes that confined the ionic motion to one dimension, and long time steps were used to save in computation time. Nevertheless, in their BD simulations of a potassium-selective channel, Bek and Jakobsson (1994) were able to capture some of the salient features of this channel. Among these are the magnitude of conductance, the shape of the current-voltage relationship, and the saturation of currents with an increasing ionic concentration.

Here we report the results of a study aimed at elucidating the permeation of ions through a model ion channel. The shape of the channel is approximated as a toroid, and cylindrical reservoirs containing 25 chloride and 25 sodium ions are placed at each end of the channel. The choice of a toroidal boundary stems from the desire to carry out the BD simulations in a realistic vestibular geometry, with all of the electric forces acting on the ions properly included. A numerical solution of Poisson's equation for an arbitrary boundary is possible, but in a BD simulation it would be computationally prohibitive. Thus one requires a coordinate system in which Poisson's equation can be solved analytically. Toroidal coordinates are the only system satisfying this dual requirement of vestibular geometry and analytical solutions (Kuyucak et al., 1997). The importance of vestibules in a channel is amply demonstrated in our simulations. For example, the repulsive dielectric force arising from the induced surface charges is strong enough to expel an ion when it is deliberately placed inside the vestibule, and a driving force provided by a membrane potential of 100 mV is not sufficient to overcome this force. Only when dipoles with a total moment of 100 × 10³⁰ Cm are placed at each end of the transmembrane segment in a favorable direction can ions drift across the channel under such a driving force. If the strengths of dipoles are increased, the average drift velocity of an ion decreases, rather than increases, owing to the resulting potential well in which the ion becomes trapped.

	METHODS

Top Abstract Introduction Methods Results Discussion References

Water-protein boundary and reservoirs

A toroidal channel was generated by rotating around the z axis two identical circles with radii of 40 Å and centers located on the x axis at 44 and 44 Å (see Fig. 1). Thus the narrowest section of the channel formed in this way had a radius of 4 Å, from which the vestibules extended to 40 Å. On each side of the vestibule, a cylindrical reservoir with radius 30 Å and length 90 Å was placed. Each reservoir contained 25 sodium and 25 chloride ions, corresponding to an ionic concentration of 150 mM. Simulations under various conditions, each lasting either 30,000 or 50,000 time steps, were repeated many times (usually nine times). For successive trials, the positions of the ions in the last time step were used as the initial starting positions of the following trial, except that a test particle was taken out and placed at a desired position.

View larger version (13K): [in this window] [in a new window]   FIGURE 1   The simulation system. A model channel with two reservoirs was generated by rotating the circles and rectangular boxes along the symmetry z axis (broken line) by 180°. The constricted region of the toroidal channel had a radius of 4 Å. Each cylindrical reservoir, 60 Å in diameter and 90 Å in height, contained 25 sodium and 25 chloride ions. Thus the ionic concentration in the volume composed of the channel vestibules and the reservoirs was 150 mM. The cylindrical reservoir had a glass boundary, in that an ion moving out of the boundary was reflected back into the reservoir.

The ion-ion potential was computed from Coulomb's law. To prevent two ions from coming too close to each other, a repulsive short-range potential, varying as 1/r⁹, was included (Pauling, 1942). This steeply rising potential imitates the repulsive force produced when the electron shells of two ions begin to overlap. An external electric field across the channel was created by applying a potential difference across the two reservoirs such that the strength of the field  would be 10⁷ V/m in the absence of the dielectric boundary.

In some simulations, a set of four dipoles was placed inside the protein boundary at z = 5 Å, and another set of four dipoles was placed at z = 5 Å. Their orientations were perpendicular to the central axis of the lumen (z axis). For each dipole, the negative pole, placed 2 Å inside the water-protein boundary, was separated from the positive pole by 5 Å. Thus, if 5/16 of an elementary charge is placed on each pole, then the total moments of four such dipoles is 100 × 10³⁰ Coulomb-meter, or ~30 Debyes (1 Debye = 3.33 × 10³⁰ Coulomb-meter). (Hereafter, Coulomb-meter will be abbreviated as Cm.) The same configuration of dipoles was used in all of the simulations, giving rise to an attractive potential for a cation and a repulsive force for an anion in the channel. We changed the total moment in the ring of four dipoles by varying the electric charges on the poles, and quote only their total strength below.

All of the above electric fields are significantly modified in the presence of the dielectric boundary. Especially when an ion is in the vestibule, induced surface charges can alter the local fields around the ion drastically, leading to an insurmountable potential barrier. The effect of the induced surface charges on the motion of ions is properly taken into account by solving Poisson's equation in toroidal coordinates, as described by Kuyucak et al. (1997).

The steep repulsive force at the dielectric boundary due to the image charges is usually sufficient to prevent ions from entering the channel protein. Nevertheless, fluctuations occur, and it is possible for an ion to gather enough momentum to overcome this repulsive force. To prevent ions from appearing inside the channel protein, we erected an impermeable glass boundary at the water-protein interface. Any ion colliding with this boundary was elastically scattered. A similar glass boundary was implemented for the reservoir boundaries.

Brownian dynamics

Brownian dynamics offers one of the simplest methods for following the trajectories of interacting ions in a fluid. The algorithm for BD is conceptually simple: the motion of the ith ion with mass m_i and charge q_i is governed by the Langevin equation

(1)

The first two terms on the right-hand side of Eq. 1 describe the effects of collisions with the surrounding water molecules. The first term corresponds to an average frictional force with a friction coefficient given by m_ii (1/_i is the relaxation time constant of the system). The second term, F_R(t), represents the random part of the collisions and rapidly fluctuates around a zero mean. The frictional and random forces in Eq. 1 are connected through the fluctuation-dissipation theorem (Reif, 1965), which relates the friction coefficient to the autocorrelation function of the random force,

(2)

where k and T are the Boltzmann constant and temperature in degrees Kelvin, respectively. Here and throughout, angular brackets denote ensemble averages. Finally, _i in Eq. 1 denotes the total electric field at the position of the ion arising from 1) other ions, 2) fixed charges in the protein, 3) membrane potential, and 4) induced surface charges on the water-protein boundary. This term in Eq. 1 is computed by solving Poisson's equation. Note that in three dimensions, Eq. 1 has to be solved for each Cartesian component (x, y, z) of the velocity.

The BD algorithm we used for solving the Langevin equation is that devised by van Gunsteren and Berendsen (1982). The steps involved in computing the positions and velocities of ions at discrete time steps of t are as follows. First, the magnitude of the electric force F(t) = q_ii acting on an ion at time t_n, and its derivative, [F(t_n)  F(t_n  1)]/t, are computed from the potential V obtained by solving Poisson's equation. Second, a net stochastic force impinging on an ion over the time period of t is computed from a sampled value of F_R(t). Finally, by substituting F(t_n), its derivative, and F_R(t) into the solutions of the Langevin equation (equations 2.6 and 2.22 of van Gunsteren and Berendsen, 1982), the position of each ion at time t_n + t and its velocity at time t_n are obtained. This process is repeated for all ions in the system for a desired number of simulation steps.

An advantage of the above algorithm over some earlier ones is that one is not limited by the condition t  1/. For the Na ions, this condition would require t  10 fs, which would severely limit the applicability of BD to ion channels. With the algorithm devised by van Gunsteren and Berendsen (1982), only two factors limit the choice of t. The average distance a particle traverses in each time step must be small compared to the dimensions of the system. Furthermore, the time derivative of the electric forces must be small relative to the absolute magnitude of the force. In a preliminary series of simulations, we have systematically increased the time step, t, from 25 fs to 1.6 ps and examined the motion of the test particle. As t increased beyond 100 fs, the trajectory of the test particle began to deviate systematically from that obtained with a short time step. In the current simulations, we have used t = 50 fs. By using the average thermal velocity for the sodium ions (see Fig. 2 B), the average displacement of ions in time t is found to be ~0.25 Å. (The radius of the narrow segment of the toroidal channel, in contrast, is 4 Å.) The change in the electric field during this time step is calculated to be a few percent at most.

View larger version (14K): [in this window] [in a new window]   FIGURE 2   Test of the Brownian dynamics algorithm. (A) The mean-square displacements, x2, of chloride ions () and sodium ions () are plotted against time. The solid lines superimposed on the graph are derived from the relation x2 = 2Dt, where the diffusion coefficients D for sodium and chloride ions are given in the Methods section. The mean square displacements determined from simulated data deviated systematically from the predicted values. The reason for these expected discrepancies between the theoretical and measured values is given in the text. (B) The Maxwellian velocities distributions, F(v)dv, of sodium () and chloride () ions are normalized for the volume of phase space and plotted. The solid lines superimposed on the measured distributions are calculated from Eq. 4. (C) The velocity autocorrelation functions, v(0)v(s), for sodium () and chloride () ions decayed exponentially. The continuous lines are calculated from Eq. 5.

The BD program was written in FORTRAN, vectorized, and executed on a supercomputer (Fujitsu VPP-300). The number of floating point operations required for each time step was on the order of 10⁹, and the number of time steps we computed was either 30,000 or 50,000. With 90% vectorization of the code, the CPU time of a supercomputer needed to complete a simulation period of 2.5 ns (50,000 time steps) was ~6 h.

The following physical constants were employed in our calculations:

Dielectric constants: _water = 78.54, _prot = 2

Masses: m_Na = 3.8 × 10²⁶ kg, m_Cl = 5.9 × 10²⁶ kg

Diffusion coefficients: D_Na = 1.33 × 10⁹ m² s¹, D_Cl = 2.03 × 10⁹ m² s¹

Relaxation time constants, ¹: _Na = 8.1 × 10¹³ s¹, _Cl = 3.4 × 10¹³ s¹

Ion radii: r_Na = 0.95 Å, r_Cl = 1.81 Å

Room temperature: T_r = 298 K

Boltzmann constant: k = 1.38 × 10²³ JK¹

	RESULTS

Top Abstract Introduction Methods Results Discussion References

Validation of the Brownian dynamics algorithm

To confirm that the behavior of the interacting ions deduced from simulations accords with the physical reality, we examined the mean square displacement, x²; the velocity distribution; and the velocity autocorrelation function, v(0)v(s). These simulations were carried out on the toroidal channel with two reservoirs, as described in the Methods section, except that the dielectric constant of the protein was set equal to that of water. Thus ions are scattered elastically from the boundary, but otherwise there were no forces acting on them because of the change in the dielectric constants.

Theoretically, the mean square displacement x² should obey the relation

(3)

In Fig. 2 A, the mean square displacement obtained for sodium (filled circles) and chloride (open circles) from one simulation lasting 500,000 time steps (25 ns of real time) is plotted against time. These are compared with the predicted slopes obtained from Eq. 3 (solid lines in Fig. 2 A). The simulation results are ~7% lower than the predicted values for a bulk solution, which is probably due to ions scattering back from the boundary, retarding their free diffusion.

Fig. 2 B shows the velocity distributions of sodium and chloride ions in the system. From the equipartition theorem, the equilibrium distribution of the velocity should be Maxwellian of the form (Reif, 1965)

(4)

where n is the number density of ions and F(v)dv is the mean number of ions per unit volume with a speed in the range between v and v + dv. The velocity distributions obtained from the BD simulations are seen to match closely those computed from Eq. 4 (solid lines).

The velocity autocorrelation function is of the form (Reif, 1965)

(5)

Thus, regardless of the initial velocity, the successive velocities will be correlated over a time interval on the order of 1/_i, the relaxation time constant of the system. To verify Eq. 5, we determined the velocities of the ions every 0.5 fs for 50,000 time steps. The autocorrelation functions shown in Fig. 2 C decayed exponentially, as predicted from Eq. 5 (solid lines). Again, filled and open circles represent the calculated values of sodium and chloride ions, respectively.

From a number of these preliminary simulations, we conclude that the BD algorithm faithfully characterizes the motion of interacting ions in a fluid confined by a glass boundary.

Repulsive dielectric force

When deliberately placed inside the toroidal channel, ions were rapidly expelled to the reservoir within 1-2 ns. A sodium ion in the system, the test particle, was placed inside the channel, on the central axis of the pore at z = 5 Å, as indicated in the inset of Fig. 3. The initial locations of the remaining 99 ions were randomly assigned. The position of the ion placed in the channel, as well as the remaining ones in the system, was calculated at each discrete time step of 50 fs for 50,000 time steps. Thus the total simulation time for one such trial was 2.5 ns. These computational steps were repeated five times. For each trial, the positions of all ions except the test particle, at the last time step in the preceding trial, were used as the initial starting positions of the subsequent trial. The same series of simulations was repeated for a chloride ion. A representative simulation from each set is shown in Fig. 3 A, where the trajectories of the sodium and chloride ions placed in the channel are plotted against time. Each of consecutive 200 points plotted in Fig. 3 A represents an average of 100 time steps or 5 ps. The ensemble averages of five trials each for the sodium and chloride ions are illustrated in Fig. 3 B. Here the positions of the ions during the entire simulation period of 2.5 ns are plotted.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Ejection of ions from the vestibule. The filled circle in the inset indicates the initial position of the input ion at the beginning of each set of simulations. Each simulation lasted 50,000 time steps. For clarity, the reservoirs and the remaining 99 ions are not shown. (A) The graphs show the positions of a sodium ion () and a chloride () ion during the first nanosecond. Each point represents the average of 100 consecutive time steps or 5 ps. From the starting position at z = 5 Å, both anion and cation rapidly moved out of the narrow segment of the channel toward the upper reservoir. (B) The mean positions at each consecutive 5 ps are averaged over five trials and then plotted against time.

The speed with which the ionic species were ejected was different. Owing to their higher mobility, chloride ions were consistently expelled faster than sodium ions. Once they were ejected from the channel, these or any other ions in the upper reservoir rarely, if ever, drifted back inside of the channel. These simulations demonstrate that the repulsive dielectric force arising from induced surface charges on the protein wall renders the channel vestibule, for the most part, devoid of ions.

The magnitude of the repulsive force presented to an ion by the dielectric boundary is sufficiently large that a driving force provided by a transmembrane potential of 100 mV cannot counteract it. This conclusion is based on the following series of simulations. After placing a sodium ion at z = 20 Å (see the inset of Fig. 4), an electric field of 10⁷ V/m was applied across the channel in the positive z direction. In the absence of any dielectric force opposing it, the ion would have drifted across the channel within 2 ns. In Fig. 4 A, three trajectories of a sodium ion are illustrated. Again, each point represents the average of 100 consecutive time steps, or 5 ps. In seven of nine trials, the ion placed at z = 20 Å either was ejected from the channel and entered the lower reservoir (filled circles in Fig. 4 A) or remained in the channel vestibule (dotted line). In the remaining two trials, the ion was able to penetrate into the other side of the channel (open circles). The average of all nine trials is shown in Fig. 4 B. The driving force provided by the applied electric field was effectively opposed by the repulsive dielectric force.

View larger version (13K): [in this window] [in a new window]   FIGURE 4   Repulsive dielectric force and applied electric field. After the test sodium ion was placed at z = 20 Å, an electric field of 107 V/m was applied in the positive z direction. Shown in A and B are, respectively, three examples of the trajectories of the test particle and the average trajectory of nine trials.

Dipoles in the transmembrane segment

The previous simulations indicate that, unless the potential barrier arising from the induced surface charges is counteracted by fixed charges or dipoles in the channel protein, ions in the majority of trials are not able to permeate the channel. In the following series of simulations, we explored the effects of such residual charges on ion permeation by placing a ring of four dipoles at each end of the transmembrane segment, as described in the Methods section (see the inset of Fig. 5). The total strength of the dipoles on each side was 100 × 10³⁰ Cm. This configuration of dipoles would be attractive to cations, but would create an additional repulsive force for anions.

View larger version (15K): [in this window] [in a new window]   FIGURE 5   Cancellation of repulsive force by dipoles. The test particle, placed at z = 5 Å, was released, and its positions during the subsequent 2.5 ns were plotted. The procedures of simulations were identical to those in Fig. 3, except that four dipoles with a total moment of 100 × 1030 Cm were placed on each side of the midline. The dipoles were placed in the channel protein such that their negative poles would face the lumen. (A) The trajectory of a chloride ion () obtained from one trial shows that it was rapidly ejected from the channel. In contrast, a sodium ion () remained in the constricted segment of the channel. (B) The trajectories for the test ions, open circles for chloride and filled circles for sodium, are obtained by averaging five successive trials.

Fig. 5 demonstrates the effects of dipoles on the motion of sodium and chloride ions that were deliberately placed inside the channel. A chloride ion placed at z = 5 Å along the central axis was expelled from the channel vestibule. In contrast, a sodium ion, when placed at the same position, remained in the vicinity of the dipoles, where charges of opposite polarity were providing an attractive potential for the ion. Sample trajectories of a chloride ion (open circles) and a sodium ion (filled circles) are shown in Fig. 5 A. The trajectories shown in Fig. 5 B are the averages of five such trials. Averaging reduced the fluctuations, making the main trend clearer. Thus a chloride ion was expelled from the channel rather quickly (in less than a nanosecond). For a sodium ion, on the other hand, the dipoles have eliminated the repulsive dielectric force that would have impinged on it, and hence it diffused freely in the channel.

Permeation of ions through the channel

The channel that had been impermeable to both cations and anions (see Fig. 4) becomes a cation-selective channel once dipoles are placed in the transmembrane segment such that their negative poles would face the channel lumen. In Fig. 6 A, we illustrate the mean positions of a sodium ion during successive 5-ps steps as it moves across the channel under the influence of an applied electric field of 10⁷ V/m. The three trajectories illustrated in the figure were obtained from the first three consecutive trials. The strength of the dipoles on each side of the transmembrane segment was 100 × 10³⁰ Cm. In all nine such attempts, a sodium ion, when released from z = 20 Å, successfully navigated the channel. The trajectory averaged over nine trials is shown in Fig. 6 B.

View larger version (14K): [in this window] [in a new window]   FIGURE 6   Permeation of sodium ions across the channel. A sodium ion was placed at z = 20 Å, as indicated in the inset, and then a potential difference of 100 mV was applied across the channel. Dipoles on each side of the transmembrane segment had a total moment of 100 × 1030 Cm. (A) Three single trajectories show that the test ions successfully traversed the channel under the influence of the driving force. (B) The trajectory illustrated is the average of nine trials. The mean drift velocity of the particles was 2.6 m s1.

The sodium ion in these simulations had an average drift velocity of 2.6 m/s in the channel, which is ~5 times faster than that of a sodium ion in bulk electrolyte solutions moving under the influence of the same applied potential gradient. In part this is because the electric field in the constricted segment of the channel, far from being constant, is enormously enhanced by the dielectric boundary and the presence of dipoles (see Kuyucak et al., 1997). Also noteworthy in Fig. 6 is the observation that the ion was not detained at the two regions where the rings of dipoles were located, indicating that the potential well created by them was not deep enough to trap the ion in it for a prolonged period of time.

The drift velocity across the channel did not increase with an increasing strength of the dipoles in the transmembrane segment. The same series of nine simulations was carried out as in Fig. 6, except that the dipole strengths were doubled and then tripled. Two typical trajectories, one obtained with a channel with a dipole moment of 200 × 10³⁰ Cm at each end of the transmembrane segment (filled circles), and the other with a dipole moment of 300 × 10³⁰ Cm (open circles), are illustrated in Fig. 7 A. The averages of nine such trials are plotted against time in Fig. 7 B, in which the trajectories shown in Fig. 4 B (no dipoles) and Fig. 6 B are reproduced for comparison. The average drift velocity of a sodium ion across the channel was reduced when the dipole strength was increased. The ion was detained for a prolonged period of time in the constricted segment of the channel. When the dipole strength in the channel was increased to 300 × 10³⁰ Cm, the ion was unable to escape from the potential well created by them during the simulation period of 2.5 ns in six of nine attempts. The conductance of the channel, however, is unlikely to be reduced by the presence of the well, because another ion of the same polarity approaching the region of the potential well will provide an additional driving force for the first ion detained near the charge residues.

View larger version (15K): [in this window] [in a new window]   FIGURE 7   Effects of increasing dipole strength on ion permeation. The trajectories of the test ions were simulated under the same conditions as in Fig. 6, except that the dipole strengths were increased first to 200 × 1030 Cm, and then to 300 × 1030 Cm. (A) Two single trajectories, obtained with two different strengths of dipoles. (B) The number accompanying each trajectory, which is the average of nine trials, represents the dipole strength (× 1030 Cm). For comparison, two average curves, indicated by 0 and 100, are reproduced from Figs. 4 and 6.

Trajectory of ions

The path taken by an ion as it journeys across the channel is predominantly along the central axis. Fig. 8 shows snapshots of the positions of a sodium ion as it traversed the channel under the influence of an electric field of 10⁷ V/m and in the presence of dipoles of strength 100 × 10³⁰ Cm (these are the same simulations mentioned in Fig. 6). Each dot in the figure represents the location of the ion on the z and x axes averaged over 100 time steps. In the first example (Fig. 8 A), the ion was slightly deflected toward the dipoles as it moved across the transmembrane segment. In the second and third examples (Fig. 8, B and C), the ions spent longer periods of time in the vicinity of the rings of dipoles, as indicated by the densities of dots. We have not expected, nor have we found, that an ion would bind to binding sites on the channel protein. These findings are in accord with those reported by Bek and Jakobsson (1994).

View larger version (16K): [in this window] [in a new window]   FIGURE 8   The path of ions. The mean positions of the test ions, as they navigate across the channel, are indicated as dots. Each dot represents the average of 100 consecutive positions of the cation as it drifts across the pore. The data were derived from those illustrated in Fig. 6. Three different examples are shown.

Shielding by mobile charges

The following series of simulations are aimed at elucidating the extent to which mobile charges are shielding permanent and induced charges on the surface of the channel protein. We placed two rings of dipoles, each with a moment of 200 × 10³⁰ Cm, at their usual positions, and then placed a chloride ion on the central axis at z = 5 Å. The speed with which such an ion was expelled from the channel was compared under three different conditions. The results of two series of simulations are shown in Fig. 9 A. In the first series of six trials (open circles), the reservoirs contained only one sodium ion. The mean positions of the test ion at successive 5-ps steps are compared with the average of six trials obtained with 50 sodium ions and 49 chloride ions in the system (solid line). That the two trajectories are broadly similar suggests that when an ion is moving out of the channel, other ions in the vicinity do not have sufficient time to rearrange and provide any discernible shielding. In the case of no dipoles in the channel (cf. Fig. 3 B), one obtains a trajectory that is very similar to those in Fig. 9 A, and the above conclusion applies.

View larger version (15K): [in this window] [in a new window]   FIGURE 9   Shielding of repulsive dielectric force by counterions. Trajectories of the test ions under various conditions were compared to assess the extent to which counterions can shield repulsive dielectric force. (A) The trajectory of a chloride ion, placed at z = 5 Å, was determined first in the presence of the remaining 99 ions in the upper and lower reservoirs (). The same series of simulations was carried out with only one sodium ion in the upper reservoir (). Each curve is the average of six trials. The strength of the dipoles on each side of the midline was 200 × 1030 Cm. (B) The test chloride ion, placed at z = 5 Å, was held at that position for 2.5 ns or 50,000 time steps, while all of the remaining ions were allowed to move. After this equilibriation period of 2.5 ns, the test ion was released, and its trajectory during the next 1.5 ns was determined (). The trajectory of the test ion, released after an equilibriation period of 10 ns, is shown by filled circles. Each curve represents the average of nine trials. The strength of the dipoles placed on each side of the midline was the same as in A. (C) A sodium ion was placed at z = 10 Å, and an electric potential of 200 mV was applied across the channel. The trajectory of the test ion, averaged over eight trials, is shown (). The trajectory shown in open cirlces, also averaged over eight trials, was obtained from the same series of simulations, except that the test particle was held at the position z = 10 Å for 2.5 ns before a potential difference was applied and it was allowed to move. The strength of dipoles on each side of the midline was 100 × 1030 Cm.

In the next two series of simulations, we address the question of equilibration. A chloride ion was placed on the central axis at z = 5 Å, and its motion was constrained, while the remaining 99 ions were allowed to move. In two series of nine trials each, this system was allowed to equilibriate for 50,000 time steps (2.5 ns) and 200,000 time steps (10 ns). Such a contrived situation was not intended to reflect a real physical situation occurring in a channel (a chloride ion would have been rapidly ejected if its motion were not constrained), but rather to see the maximum possible effect of shielding when all of the forces among the ions and the boundary were properly taken into account. The results of these simulations are shown in Fig. 9 B, together with the no-equilibriation case from Fig. 9 A. There is negligible difference in the motion of the test ion between the two equilibriation periods, which indicates that the period was sufficiently long. Compared to the no-equilibriation case, the test ion moves relatively slowly, which points to some shielding effect. But it is not large enough to change the dynamics of the ion; the ion is still ejected, albeit somewhat more slowly. To make these statements on shielding more quantitative, we compare the average net force the test ion experienced during the first 0.25 ns (or 5000 time steps). When the reservoirs contained only one counterion, this force was 10.5 pN. The equivalent force experienced by a chloride ion when the reservoirs contained 150 mM sodium and chloride ions and the ion was released after the equilibriation period of 10 ns was 8.1 pN. Thus the repulsive force was reduced at most by 23% when the test ion was forcibly held for a prolonged period of time, allowing counterions to come near.

Finally, to mimic a realistic situation, two rings of dipoles, each with a moment of 100 × 10³⁰ Cm, were placed in the channel as before, and a sodium ion was held at z = 10 Å. The dipoles whose negative poles faced the channel lumen would have canceled the repulsive dielectric force impinging on a sodium ion, and its presence in the channel would not have created any energy well to attract counter-ions. After applying a potential difference of 200 mV, the positions of the test particle in the following 1.5 ns were determined. The trajectory of the test particle, averaged over 8 trials, is shown Fig. 9 C (filled circles). The same series of simulations was carried out as before, except that the test particle was held at z = 10 Å for 2.5 ns before the potential difference was applied across the channel. During this equilibriation period, all of the remaining 99 ions were allowed to execute Brownian motions. The trajectory of the test ion after equilibriation, averaged over eight trials (shown as open circles), is compared with the trajectory obtained without equilibriation. The two curves do not differ appreciably.

The results of these simulations demonstrate that shielding does not play a significant role in channel dynamics. Even when its effects were artificially maximized (see Fig. 9 B), the magnitude of dielectric force impinging on an ion was reduced only slightly.

	DISCUSSION

Top Abstract Introduction Methods Results Discussion References

The Brownian dynamics algorithm devised by van Gunsteren and Berendsen (van Gunsteren and Berendsen, 1982; van Gunsteren et al., 1981) is highly suitable for studying the motions of interacting ions in a fluid. In this algorithm, the stochastic force featured in the Langevin equation is treated such that the range of time steps t one is permitted to choose is not restricted by the relaxation time constant of the system. Moreover, the integration of the electric force is performed in a way similar to that of the molecular dynamics algorithm of Verlet. Indeed, the algorithm reduces to the Verlet algorithm when the friction coefficient is reduced to zero, and moreover, it traces diffusion of interacting ions when the time step t is made large relative to the relaxation time constant. Our preliminary tests of this algorithm in a simple periodic boundary and a toroidal glass boundary showed that the computational steps capture many of the important features of the motions of interacting ions in a fluid. Among these are the mean square displacement (Fig. 2 A), the velocity distribution (Fig. 2 B), and the characteristic relaxation time constant (Fig. 2 C). Moreover, the temperature of the system remained stable over the simulation period, and the ionic concentration in a localized region, revealed by the time average of the probability of occupancy, remained constant at 150 mM.

With this algorithm we were able to unravel several salient features of ion permeation across a model biological channel. First, the repulsive dielectric force presented to an ion by the vestibular wall is large (Fig. 3). The driving force provided by a membrane potential of 100 mV is insufficient to counteract this force (Fig. 4). Second, fixed negative charge residues (or dipoles) near the transmembrane segment render the channel permeable to cations, while enhancing the potential barrier seen by anions (Fig. 6). Third, the presence of an excessive number of dipoles impedes the passage of an ion across the channel (Fig. 7). This is because an ion traversing the channel is detained in a deep potential well created by the charge residues and is unable to climb out of the well until it acquires either a large kinetic energy or another ion of the same polarity approaches the detained ion. Finally, we show that the trajectory taken by an ion as it navigates the channel is primarily along the central axis of the pore (Fig. 8). This finding is in accord with the inference drawn by Bek and Jakobsson (1994) from the results of their BD simulations in the diffusive regime.

The behavior of an ion in the channel can be inferred qualitatively by constructing a potential profile obtained by placing an ion at various locations (Hoyles et al., 1996; Kuyucak et al., 1997). Such potential profiles, however, do not enable us to make accurate predictions about the motions of interacting ions in the ensemble at a microscopic time scale. Ions and water molecules continuously collide with each other under thermal motion. In BD calculations, all of these molecular impacts impinging on an ion are lumped together into a stochastic force. The magnitude of this force is large compared to the electric force experienced by an ion at any given time, as is the thermal velocity compared to the drift velocity. Moreover, a potential profile merely tells us the net electric force experienced by an ion in the absence of any other ions in the vicinity. To determine the position of an ion during a discrete time step or to trace its trajectory in a fixed time interval, both stochastic force acting on an ion at a given time and the Coulomb forces arising from all of the other charges in the system must be included as part of the total force influencing its motion. Many of the features we have uncovered from the BD simulations of the permeation of ions are consistent with the intuitive picture obtained from the potential profiles as illustrated in Kuyucak et al. (1997).

The extent to which the ionic atmosphere created by mobile charges can shield fixed charges in the protein and induced charges on its surface has been addressed by other investigators (Jordan, 1987; Jordan et al., 1989; Chen et al., 1997). The result we report here (Fig. 9) is the first study to assess in detail the contributions made by mobile ions on the net force experienced by an ion in a model channel that approximates the shape of a real biological channel. Ions in our simulations execute thermal motions in a three-dimensional space, and the average distance traversed by them in each time step of 50 fs is ~0.25 Å. The small time step used relative to the mean free path ensured that ion-ion and ion-protein distances closely mirror the real situation in electrolyte solutions. Our analyses reveal that counterions have a negligible influence on shielding of the dielectric and dipole forces on an ion in the model channel.

The role of the vestibule in the permeation of ions across a membrane channel can be summarized as follows. A large vestibule renders the channel impermeable to ions. Fixed charges or dipoles near the transmembrane segment of the protein eliminate the potential barrier seen by one ionic species while enhancing the barrier seen by the other. A channel thus becomes either cation- or anion-selective, and the high potential barrier ensures that counterions do not leak across the opposite direction under the influence of an applied electric field. It is important to realize that a vestibule of the type imaged by Toyoshima and Unwin (1988) contains ~500 water molecules and, on average, one cation at 150 mM, provided the dielectric force is completely canceled by fixed negative charges. A second cation approaching the vestibule will induce surface charges on the protein wall, which will tend to repulse the ion. In addition, there will be a Coulomb repulsive force between two ions; these effects combined effectively preclude the entry of a second or third ion into the vestibule. The mean ion-ion distance in the vestibule under various conditions, together with its implication for the conductance-concentration relationship of single channels, will be detailed in a subsequent report.

We have emphasized previously (Hoyles et al., 1996; Kuyucak et al., 1997) and reiterate here that inferences made from macroscopic electrostatics are valid only in the regions that are large compared to the diameters of water and ion molecules. In the narrow constricted region of the channel, whose radius may be only ~4 Å (nearly equal to that of a potassium ion with its first hydration shell), the representation of the dielectric as a continuous medium is a poor approximation. Here polar groups on the wall will interact with water molecules surrounding an ion as it attempts to traverse the pore, and the ability of resident water molecules to align with an electric field will be severely restricted. It is most likely that this interaction between the dipoles on the protein wall and the primary or secondary hydration shells of an ion is what determines the ionic selectivity of anion or cation channels. Studies of the temperature dependence of the conductivity of ion channels reveal that ions traversing the membrane encounter an energy barrier of ~3kT_r, which suppresses the current by an order of magnitude from that predicted from the geometrical considerations alone (Kuyucak and Chung, 1994; Chung and Kuyucak, 1995; Milburn et al., 1995). This barrier may well arise from the dynamic interactions taking place in the constricted channel segment between the protein wall and the hydration shells. Elucidation of the permeation process taking place in this region requires molecular dynamics calculations such as those reported for the gramicidin channel (e.g., Mackay et al., 1984; Roux and Karplus, 1991; Chiu et al., 1989).

This work was supported by grants from the Australian Research Council and the National Health and Medical Research Council of Australia. The calculations upon which this work is based were carried on the Fujitsu VPP-300 of the ANU Supercomputer Facility.