| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Biophys J, August 1998, p. 793-809, Vol. 75, No. 2
*Protein Dynamics Unit, Brownian dynamics simulations have been carried out to
study ionic currents flowing across a model membrane channel under various conditions. The model channel we use has a cylindrical transmembrane segment that is joined to a catenary vestibule at each
side. Two cylindrical reservoirs connected to the channel contain a
fixed number of sodium and chloride ions. Under a driving force of 100 mV, the channel is virtually impermeable to sodium ions, owing to the
repulsive dielectric force presented to ions by the vestibular wall.
When two rings of dipoles, with their negative poles facing the pore
lumen, are placed just above and below the constricted channel segment,
sodium ions cross the channel. The conductance increases with
increasing dipole strength and reaches its maximum rapidly; a further
increase in dipole strength does not increase the channel conductance
further. When only those ions that acquire a kinetic energy large
enough to surmount a barrier are allowed to enter the narrow
transmembrane segment, the channel conductance decreases monotonically
with the barrier height. This barrier represents those interactions
between an ion, water molecules, and the protein wall in the
transmembrane segment that are not treated explicitly in the
simulation. The conductance obtained from simulations closely matches
that obtained from ACh channels when a step potential barrier of 2-3
kTr is placed at the channel neck. The
current-voltage relationship obtained with symmetrical solutions is
ohmic in the absence of a barrier. The current-voltage curve becomes
nonlinear when the 3 kTr barrier is in
place. With asymmetrical solutions, the relationship approximates the
Goldman equation, with the reversal potential close to that predicted
by the Nernst equation. The conductance first increases linearly with
concentration and then begins to rise at a slower rate with higher
ionic concentration. We discuss the implications of these findings for
the transport of ions across the membrane and the structure of ion
channels.
Theoretical studies of the biological ion channel
have been hampered by a lack of detailed structural knowledge. The
exact shape of any biological channel, either ligand-gated or
voltage-activated, is unknown, as are the positions, densities, and
types of dipoles and charge moieties on the protein wall. These details
are needed to compute the intermolecular potential operating between
water molecules, ions, and the protein wall, which is the essential ingredient for microscopic studies of channels using molecular dynamics. Even if such information were available, it would not be
feasible at present to carry out molecular dynamics calculations for
all of the water molecules and ions in a biological channel and its
surroundings for a period long enough to deduce any of its
macroscopically observable properties. For these and other reasons, it
has not been possible to compare the results obtained from microscopic
models with the real data obtained from patch-clamp recordings. There
is a need to develop models that can relate the structural parameters
of channels to experimental measurements and to build a theoretical
framework that interlaces all of the disparate sets of observations
into a connected whole. The theoretical model described in this
paper was produced in the hope of furthering this aim.
It is possible to make the computations tractable by making several
simplifications of and idealizations about the channel and electrolyte
solutions, and to examine the magnitude of currents flowing through the
model conduit under various conditions, using Brownian dynamics
simulations (Cooper et al., 1985 To deduce the conductance of a model channel with Brownian dynamics
simulations, the number of ions placed in the reservoirs that mimic the
extracellular and intracellular media must be relatively large, and the
simulation period needs to be sufficiently long that reliable
statistics for currents flowing across the channel can be collected.
With this requirement in mind, we have devised a method of reducing the
amount of computational effort involved in simulating a system of
charged particles interacting with a protein wall. Instead of computing
the force acting on an ion at each position and at each time step, we
precalculated the values of the electric potential and field for a grid
of positions and stored them in a set of lookup tables. Using a
multidimensional interpolation algorithm (Press et al., 1989 Brownian dynamics simulations do not adequately capture the transport
process of ions in the narrow, constricted channel region. This is
because water is treated as continuum, ions are idealized as point
particles, and the protein wall is represented as a structureless, rigid, and smooth dielectric surface. In this narrow cylindrical region, polar groups on the protein wall are likely to interact with
the primary or secondary hydration shell of an ion as it drifts across
the conduit, possibly replacing several water molecules in the shell
and forming temporary hydrogen bonds with the ion. Elucidation of the
permeation processes taking place in the transmembrane segment will
require molecular dynamics calculations, such as the ones carried out
for the gramicidin pore (Roux and Karplus, 1991a Here we describe a model channel and report the results of Brownian
dynamics calculations aimed at elucidating the permeation of ions
through it. The shape of the channel is made approximately the same as
that of the ACh channel (Toyoshima and Unwin, 1988 The channel model
There are two major impediments to formulating a microscopic model
of membrane ion channels and simulating molecular trajectories of
idealized particles interacting with the protein boundary. First, the
description of the protein wall forming a narrow cavity is relatively
incomplete. For example, the precise shape, and the number, magnitude,
and location of dipoles or charge moieties on the protein wall of
biological ion channels are unknown. Second, even with the most
advanced computer currently available, it is not feasible to simulate
motion of all water molecules and ions in a channel and surroundings
for a period long enough to measure conductance. To study the
permeation of ions through a biological channel, we need to devise an
idealized model channel and then make several simplifying assumptions
about the intermolecular potential that operates between particles in
the assembly. We therefore make the following simplifications and
idealizations about our model channel to enable us to simulate current
flow across the channel under various conditions.
Shape of the channel
ABSTRACT
Top
Abstract
Introduction
Methods
Results
Discussion
Concluding Remarks
References
INTRODUCTION
Top
Abstract
Introduction
Methods
Results
Discussion
Concluding Remarks
References
). Brownian dynamics provides one of
the simplest methods for following the trajectories of idealized
particles in a fluid interacting with a dielectric boundary. Here the
water is treated as a continuum, and the motions of individual ions are
assumed to be governed by electrostatic forces emanating from other
ions, fixed charges in the proteins, the applied electric field and the
dielectric boundary. The effects of solvation and the structure of
water are taken into account by frictional and random forces acting on
ions. Brownian dynamics simulations have already been fruitfully utilized to investigate the movement of ions across a model cylindrical tube (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and
Jakobsson, 1994) and a toroidal channel (Li et al., 1998). With these
simulations, it was possible to capture some of the salient features of
biological ion channels and reveal the importance of the vestibules in
influencing the permeability properties of ions through the pore.
), the
field and potential experienced by an ion at any position could be
deduced from the information stored in the lookup tables. Using this
method, we were able to study ionic currents flowing across a realistic
model channel under various conditions.
). In the absence of a
detailed knowledge of the location and types of polar groups lining the
channel wall and the precise geometry of the transmembrane segment, we
have represented the intermolecular interactions taking place between
ions, water molecules, and the protein wall in this region as a step
potential barrier of variable height. The barrier is constructed such
that it rises gently to the desired height in 1 Å, and its first and
second derivative are zero at the end points. Thus the energy must be paid to enter the neck and is returned when the ion exits.
), and cylindrical
reservoirs containing sodium and chloride ions are placed at each end
of the channel. We show that by placing an appropriate strength of
dipoles in the channel wall and erecting a small energy barrier at the
entrance of the narrow transmembrane segment, we can replicate some of
the macroscopically observable properties of biological ion channels.
Among these are the channel conductance, an ohmic current-voltage
relationship, inward rectification as predicted by the Goldman
equation, and the conductance-concentration relationship.
METHODS
Top
Abstract
Introduction
Methods
Results
Discussion
Concluding Remarks
References
), was
generated by a hyperbolic cosine function, z = a cosh x/a, where a = 4.87 Å. The radius of the entrance of the vestibule was fixed at 13 Å. Two identical vestibules were connected with a cylindrical
transmembrane region of length 10 Å and radius 4 Å. A cross section
of the model channel, the total interior volume of which was 2.16 × 10
26 m3, is illustrated in Fig. 1
B. We assumed, for convenience, that the two vestibules are
identical in size, although the image of the channel produced by
Toyoshima and Unwin (1988) shows that the extracellular vestibule is
larger than the intracellular vestibule.
View larger version (32K):
[in a new window]
FIGURE 1
Idealized biological ion channel. A model channel with
two catenary vestibules was generated by rotating the closed curves
outlined in A along the symmetry z axis
(horizontal line) by 180°. A transverse section of the
channel so generated is shown in B. To approximate the
shape of the acetylcholine receptor channel, vestibules at each side of
the membrane were constructed by using a hyperbolic cosine function,
z = a cosh
x/a, where a = 4.87 Å. The radius of the entrance of the vestibule was fixed at 13 Å, and
the cylindrical transmembrane segment had a radius of 4 Å. Each
cylindrical reservoir, 60 Å in diameter and 22 Å in height, contained
a fixed number of sodium and chloride ions. Unless stated otherwise,
the ionic concentration in the volume composed of the channel
vestibules and the reservoirs was 300 mM. The cylindrical reservoir had
a glass boundary, in that an ion moving out of the boundary was
reflected back into the reservoir.
Dipoles in the protein wall
To investigate how the permeation of ions across the channel is influenced by the presence and absence of dipoles in the protein wall, we placed in some simulations a set of four dipoles inside the protein boundary at z = 5 Å and another set of four dipoles 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 was placed on each pole, then the total moments of four such
dipoles would be 100 × 1030 Coulomb-meter. The same
configuration of dipoles was used in all of the simulations, giving
rise to an attractive potential for sodium ions and a repulsive
potential for chloride ions in the channel. These fixed charges
represent the charged side chains thought to form a ring around the
entrance of the constricted region (Unwin, 1989), and their nearby
counter-charges. For convenience, we adjust the amount of charge rather
than the number or positions of the charges, but in reality the side
chains of charged amino acid residues would have one electron charge
each if they are fully ionized or unprotonated at neutral pH. Polar
groups located in the transmembrane segment of the channel that may
rotate in and out form temporary hydrogen bonds with an ion navigating
across it, as found in gramicidin A pores, are not explicitly modeled in this study.
Energy barrier to penetrating the transmembrane segment
An ion in the vestibule needs to surmount an energy barrier to traverse the narrow, constricted segment of the channel. The presence of such an additional energy barrier in the gramicidin pore has been revealed by molecular dynamics calculations. Additional evidence is provided by the conductance-temperature curves measured from biological ion channels, which are always steeper than the conductivity-temperature curve of bulk electrolyte solutions. This steeper rise in the channel conductance with temperature is interpreted as being due to the presence of a dynamic barrier that ions need to overcome (Kuyucak and Chung, 1994; Chung and Kuyucak, 1995).
Intuitively, this barrier arises from the interactions between the
protein wall and the hydrated ion as the ion negotiates its way into
the narrow, cylindrical transmembrane pore. To enter the narrow
segment, the hydration shell of an ion needs to be rearranged, or some
of the water molecules in the primary or secondary hydration shell need
to be substituted with polar groups on the protein wall. To rearrange
the ion-water geometry requires an additional energy, and the ion can
surmount such a barrier only when it gains a sufficient kinetic energy.
To mimic a barrier present in the ion channel, we placed in some
simulations a potential step near the constricted segment of the
channel. The method we used for implementing such a potential barrier
in the Brownian dynamics algorithm is detailed in the Methods.
Water as a continuum
We treat the water as a continuum in which the ions move under the influence of electrostatic forces and random collisions. In the constricted region of the channel, where the radius is ~4 Å, the representation of water molecules as a continuous dielectric medium is a poor approximation. The addition of energy barriers to the model is an attempt to compensate for this difficulty. We represent the water-protein interface as a single sharp and rigid boundary between dielectrics. In reality, however, the channel wall is not made of a structureless dielectric material. Instead, its surface is likely to be lined with hydrophilic and polar side chains, which will restrict the ability of water molecules to align with the electric field. The interface could be more accurately represented as a region of intermediate dielectric constant. Because the use of a single boundary does not introduce significant error, as has been shown elsewhere (Hoyles et al., 1996), we adopt the simpler model. The polar groups and
ordered water at the interface are not explicitly included in our
model, being represented by the dielectric boundary.
Applied electric field
A potential difference across a lipid membrane is produced by a surface charge density on each side of the membrane. In microscopic terms, the surface charge density is a cloud of unpaired ions on either side of the membrane. Because these clouds are too diffuse to be explicitly included in our simulation, we apply an external electric field
of a constant strength to represent the average effect of the
ionic clouds. In the absence of any dielectric boundary, the potential
difference across a channel with the length d is /d. The presence of a dielectric boundary, however,
severely distorts the field, enhancing it in the transmembrane segment and attenuating it in the vestibule (Kuyucak et al., 1998). Thus the
precise potential difference will depend on the selected reference points at the two sides of the catenary channel. For simplicity, we
apply a field strength of 107 V m1 and refer
to it as an applied potential of 100 mV.
Brownian dynamics
The trajectories of ions drifting across the channel under the
influence of a driving force while executing random thermal motion were
followed with Brownian dynamics simulations. A detailed account of the
theory underlying Brownian dynamics calculations is given elsewhere (Li
et al., 1998). Briefly, the computational method traces the motion of
the ith ion with mass mi and charge qi with the Langevin equation:
|
(1) |
i, where 1/i is
the relaxation time constant of the system. The second term,
FR(t), represents the random part of
the collisions and rapidly fluctuates around a zero mean. The
frictional and random forces in Eq. 1, which together describe the
effects of collisions with the surrounding water molecules, are
connected through the fluctuation-dissipation theorem (Reif, 1965),
which relates the friction coefficient to the autocorrelation function of the random force,
|
(2) |
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.
Computational steps
The Brownian dynamics algorithm we used for solving the Langevin
equation is that devised by van Gunsteren and Berendsen (van Gunsteren
and Berendsen, 1982; van Gunsteren et al., 1981). That the behavior of
the interacting ions deduced from simulations using this algorithm
accords with the physical reality is detailed elsewhere (see Li et al.,
1998). To compute the positions and velocities of ions at discrete time
steps of 
t, the algorithm takes the following
computational steps:
Step 1. Compute the magnitude of the electric force
F(t) = qii acting on each ion at time
tn and calculate its derivative [F(tn) F(tn 1)]/t.
Step 2. Compute a net stochastic force impinging on each ion over the
time period of t from a sampled value of
FR(t).
Step 3. Determine the position of each ion at time
tn + t and its velocity at time
tn by substituting
F(tn), its derivative F'(tn) and
FR(t) into the solutions of the
Langevin equation (equations 2.6 and 2.22 of van Gunsteren and
Berendsen, 1982).
Step 4. Repeat the above steps for the desired number of simulation steps.
Throughout, we used the time step of t = 100 fs, except when an ion entered one of the two imaginary bands placed
around the rising and falling edges of the step energy barrier (see the following section).
Implementation of a step potential barrier in the algorithm
The use of a long time step causes a problem in implementing
potential barriers and steps in Brownian dynamics as short-range forces. In our simulations, the Brownian dynamics algorithm operates predominantly in the diffusive regime. In other words, random forces
are far more important than the velocity on the previous step in
determining the anion's new velocity and position. The algorithm
devised by van Gunsteren and Berendsen (1982) uses the factors
[exp(t)] and [1 exp(t)] to switch between kinetic and diffusive
regions. With the long time step that we use (t = 100 fs),
|
|
12 N (3 kTr over 1 Å) produces a drift velocity of 39 ms1 for sodium and a displacement of 0.04 Å in one time
step. This is only 1/7 of the average random displacement in one time
step
an ion can diffuse right through a potential barrier before the
barrier force has time to act.
To obviate these problems, we used two different time steps: a short time step of 1 fs when an ion was in the process of climbing or descending the barrier, and a long time step of 100 fs otherwise. A smooth potential barrier of height VB was erected at z = ±10 Å (5 Å from the entrance of the cylindrical segment), with the profiles at the ends as
|
(3) |
|
(4) |
z = 1 Å is the width of
the profile. The potential profile U(s) is chosen
such that it rises from zero at z = ±10.5 to
VB at z = ±9.5, and the first
and second derivatives of U(s) vanish at
z = ±9.5 and z = ±10.5. The force due
to the barrier is obtained by differentiating Eq. 3.
We included a safety distance of 0.5 Å in the potential profile. Thus,
whenever ions were in the band of z = 9 to 11 Å or z = 9 to 11 Å, we switched from long time steps to
an equivalent sequence of short time steps for those ions. Trajectories
of ions in these bands were determined by a sequence of 100 short time steps for the subsequent 100 fs. In the meantime, all of the other ions
were simulated for a single long time step in the normal way. The
long-range (electrostatic) forces on the ion were calculated at the
start of the sequence of 100 short steps and held constant thereafter.
Similarly, reflection from the boundary walls was done once at the end
of the sequence of short steps. The force from the barrier, however,
was recalculated for each short step, thus ensuring that the effect of
the barrier on the ion's motion would be accurately simulated.
Lookup tables for the electric field and potential
To reduce the amount of computational effort, we have made use of
lookup tables. The electric field and potential on a grid of positions
were precalculated, using a numerical method of solving Poisson's
equation (Hoyles et al., 1996), and then recorded in a set of tables.
By interpolating between table entries, the electric field and
potential anywhere in and in the vicinity of the channel can be quickly
determined. Technical details of this method will be described
elsewhere (Hoyles et al., 1998).
The electric potential V experienced by an ion is composed of the following four parts:
|
(5) |
experienced by an ion
in the same way:
|
(6) |
To allow rapid lookup, the precalculated results must be on an evenly
spaced grid. Because the use of a rectilinear grid would result in many
wasted points and a jagged edge near the pore boundary, we made use of
a system of generalized cylindrical coordinates for constructing the
lookup tables. The position of an ion was first converted to
generalized coordinates, and then the values of the electric potential
and field in the vicinity of the ion were extracted from the tables by
multidimensional linear interpolation, which is a simple algorithm that
generalizes easily to dimensions greater than 2 (Press et al., 1989).
We used separate tables for the self-potential
(VS,i), the external potential
(VX,i), and the image potential
(VI,ij). These are, respectively, two-, three-, and five-dimensional tables. Cylindrical symmetry has been exploited to
reduce by one the number of dimensions of the self-potential and
image-potential tables. The electric field was stored and extracted in
the same way as the potential, except that three values were required
for each point in a table, one for each Cartesian component of the
field. So while the field tables were indexed by generalized
coordinates, their contents were stored as Cartesian coordinates.
Illustrated in Fig. 2 are a graph of
potential energy (A) and that of the z component
of force (B) for a single ion moving parallel to the central
axis of a catenary channel with no dipoles. The solid lines are
calculated by using an iterative method (Hoyles et al., 1996), and the
circles by interpolating from the precalculated values stored in the
lookup tables. The spacing between points in the lookup table is 1.77 Å in the z direction, and the circles are at the midpoints
of these intervals, where the maximum interpolation error is expected
to occur. There are 37 points across the lookup table in the
r direction (perpendicular to the z axis), and
the spacing between these varies with the radius of the channel. As indicated in the inset of Fig. 2, the path of the ion is offset from
the z axis by 3 Å, and the ion passes through the midpoint of the radial interpolation intervals in the neck region.
|
The correspondence between the iterative and lookup methods evident in
Fig. 2 indicates that interpolation error is negligible for potential
energy and the z component of force in the most important
parts of the channel, namely, the center of the vestibule and the
neck region. More detailed tests indicate that the relative error
increases when an ion approaches the vestibular wall. For example, at 1 Å from the wall, the relative error in the repulsive force is ~15%.
This is not of great concern, because ions in the vestibule tend to
stay away from the water-protein boundary (Li et al., 1998). In the
constricted region, where ions are forced into closer proximity with
the wall, the error in the repulsive force 1 Å from the wall is
~5%.
Input parameters and details of simulations
Simulations under various conditions, each lasting between 500,000 and 2,000,000 time steps, were repeated many times, for five to nine trials. For the first trial, the positions of ions in the reservoirs were assigned randomly with the proviso that the minimum ion-ion distance should be 2.7 Å, or 1.5 times the radius of a chloride ion. We note that with all subsequent ion-ion distance checking, to be explained below, the minimum allowed distance, which we refer to as the "safe distance," was chosen to be 3/4 of the sum of the ionic radii. For successive trials, the positions of the ions in the last time step were used as the initial starting positions of the following trial. The current (in pA) was extrapolated from the total number of ions traversing the channel over the simulation period.
On each side of the vestibule, a cylindrical reservoir with radius 30 Å and an adjustable height was placed. A fixed number of sodium and chloride ions were placed in each reservoir, and the height of the cylindrical reservoir was adjusted to give a desired ionic concentration. As ions were forbidden to approach the wall of the reservoir within 1 Å, the effective radius of the cylindrical reservoir was 29 Å. For example, if 13 sodium and 13 chloride ions were placed in each reservoir and the desired ionic concentration was 300 mM, the height of each of the two cylindrical reservoirs was adjusted to 22 Å.
To prevent two ions from coming too close to each other, a repulsive
short-range potential, varying as 1/r9, was
included (Pauling, 1942). This steeply rising potential imitates the
repulsive force produced when the electron shells of two ions begin to
overlap. When the ionic concentration in the reservoirs was high, ions
at times were able to jump large distances and end up very close to
another ion. The forces at the next time step in such instances would
be very large, and the affected ions could leave the system. To correct
this problem, we checked ion-ion distances at each time step. If two
ions were closer than the defined safe distance, then their
trajectories were traced backward in time until such a distance was
exceeded. By performing these checks and corrections, the system was
well behaved over the simulation, even for very high concentration. Such a minor readjustment of the position of an ion was needed about
once every 100 time steps when the reservoirs and the channel contained
52 ions. The steep repulsive force at the dielectric boundary due to
the image charges was usually sufficient to prevent ions from entering
the channel protein. We ensured that no ions would appear inside the
channel protein by erecting 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.
To ensure that the desired intracellular and extracellular ion concentrations were maintained throughout the simulation, a stochastic boundary was applied. When an ion crossed the transmembrane segment, an ion of the same species was transplanted so as to maintain the original concentrations on both sides of the membrane. For example, if a sodium ion from the left-hand side of the channel crossed the narrow transmembrane segment and reached the imaginary plane at z = 10 Å, then the furthermost sodium ion in the right-hand reservoir was taken out and placed in the far left-hand side of the left reservoir. When transplanting ions, we chose a point no closer to another ion than the defined safe distance. The stochastic boundary trigger points, located at z = ±10 Å, were checked at each time step of the simulation.
The Brownian dynamics program was written in FORTRAN, vectorized, and executed on a supercomputer (Fujitsu VPP-300). The amount of vectorization varied from 67% to 92%, depending on the number of ions in the reservoirs. With 52 ions in the reservoirs, the CPU time of a supercomputer needed to complete a simulation period of 1.0 µs (10 million time steps) was ~18.7 h. The following physical constants were employed in our calculations:
Dielectric constants: 
water = 80, prot = 2
Masses: mNa = 3.8 × 1026 kg,
mCl = 5.9 × 1026 kg
Diffusion coefficients: DNa = 1.33 × 109 m2 s1,
DCl = 2.03 × 109 m2 s1
Relaxation time constants, 1: Na = 8.1 × 1013 s1, Cl = 3.4 × 1013 s1
Ion radii: rNa = 0.95 Å, rCl = 1.81 Å
Room temperature: Tr = 298 K
Boltzmann constant: k = 1.38 × 1023
J K1
Elementary charge: e = 1.60 × 1019
C
Throughout we give energy in temperature units,
kTr. We note that 1 kTr
equals 4.11 × 1021 J and 2.478 kJ/mol. Units of
dipole moments are quoted in Coulomb-meters, abbreviated hereafter as
Cm. One Debye corresponds to ~3.33 × 1030 Cm.
| |
RESULTS |
|---|
|
Top
Abstract
Introduction
Methods
Results
Discussion
Concluding Remarks
References
|
|---|
Dipoles in the channel
In the absence of any charge moieties or dipoles on the protein
wall, the potential barrier presented to an ion moving under the
influence of an applied potential of 100 mV is shown in Fig. 3 (top curve, labeled 0). The
potential energy of a sodium ion placed at a fixed position on the
z axis was calculated numerically by solving Poisson's
equation iteratively, as detailed previously (Hoyles et al., 1996). The
ion was then moved by 1 Å and the calculations were repeated. The
potential profile presented to the ion as it moved from outside
(left-hand side) to inside (right-hand side) increased slowly, peaking
at the center of the cylindrical transmembrane segment (labeled 0 Å),
and then decreased steadily as it traversed the second half of the
channel. Note that without the membrane potential, one would have
obtained a symmetrical, bell-shaped barrier with a peak height of
14.5 × 1021 J. The presence of the membrane
potential had lowered the relative height of the barrier to 3.5 × 1021 J and distorted the shape of the profile to an
asymmetrical curve.
|
Two rings of dipoles, together with an applied electric potential of
100 mV, eliminated the repulsive dielectric force. The number
accompanying each curve in Fig. 3 represents the total strength of four
dipoles (×1030 Cm) in each ring. Because there were two
rings of dipoles, one at z = 5 Å and the other at
z = +5 Å, the strength of dipoles placed on the entire
channel wall was twice the value indicated in the figure. In the
presence of dipoles, an ion traversing from outside to inside would
encounter an attractive potential throughout its trajectory. With each
stepwise increase in dipole strength, what used to be a potential
barrier became a potential well whose depth increased with the strength
of dipole moments.
The potential profiles illustrated in Fig. 3 were constructed under the assumption that there were no other ions in the system. Such profiles merely reveal that the permeation of ions across the channel would be hindered if no dipoles were present on the channel wall. It is not possible, however, to deduce the conductance of the channel from such profiles, or the effect a deep potential well will have on ions permeating the channel. To study these properties, one needs a dynamical theory.
We have used Brownian dynamics simulations to deduce the conductance of
the model channel under various conditions. In all of the subsequent
figures, unless stated otherwise, each point is the average of nine
simulations, each lasting for 500,000 time steps (50 ns), or five
simulations, each lasting for 2,000,000 time steps (200 ns). The error
bar accompanying each data point is one standard error of mean. The
error bar is not shown if it is smaller than the size of the data
point. Again, unless noted otherwise, 13 sodium and 13 chloride ions
were placed in the left-hand reservoir (representing the extracellular
space), whose volume was 5.84 × 1026
m3, and the same number of ions in the right-hand side
reservoir (the intracellular space). Thus the ionic concentration in
the reservoirs was 300 mM, which is about double the physiological concentration. This higher concentration was preferred in the simulations to obtain better statistics.
In Fig. 4, the number of ions that
traversed the channel under the driving force of 100 mV during a
simulation period of 0.45 µs (4.5 million time steps) is converted to
current in pA and plotted against the dipole strength. On the
right-hand ordinate of Fig. 4, current in pA is converted to
conductance in pS at the physiological concentration of 150 mM. Because
the current in these ranges of ionic concentrations increases almost
linearly with concentration (see later), such an extrapolation of
conductance results is justified. With no dipoles placed on the channel
wall, the number of sodium ions that traversed from outside to inside in 0.45 µs was 10, which corresponds to a current of 3.6 pA. The net
current increased rapidly with the increasing dipole strength up to
100 × 1030 Cm (filled circles in Fig.
4). The number of ions crossing the channel increased further with a
further increase in the dipole strength, but many ions were also
traversing in the opposite direction, against the direction of the
applied electric field. The current flowing from inside to outside is
indicated as open circles in Fig. 4. As a result, the net current
actually decreased as the dipole strength was increased further from
300 to 600 × 1030 Cm.
|
The strength of dipoles in biological ion channels is likely to be
sufficiently large to cancel the repulsive dielectric force presented
to the ion, but not so strong as to allow ions to move against the
applied electric gradient. If we assume that the channel shape can be
idealized as our model channel and that the dielectric constant of
water in the vestibule is 80, then the optimal strength of dipoles on
the channel wall is between 100 and 200 × 1030 Cm.
In reality, the total dipole moment present on the channel wall will be
greater than the value we quote, because the dielectric constant of
water in the vestibule must almost certainly be lower than that in bulk
water (see, for example, Gutman et al., 1992; Sansom et al. 1997).
Hereafter, we place two rings of dipoles, each with a moment of
100 × 1030 Cm, just above and below the
transmembrane segment to investigate other macroscopically observable
properties of membrane ion channels.
Energy barrier near the transmembrane segment
It is not known how large an energy barrier an ion must overcome
to cross the transmembrane segment of the channel. From the conductance-temperature relationship, it has been estimated that the
height of this barrier is ~3 kTr (Kuyucak and
Chung, 1994). Here we examine how the height of such a barrier
influences the conductance of the channel.
The number of ions crossing the channel was tabulated at a given height
of this step potential barrier during the simulation period of 1.0 µs. The current was not attenuated appreciably when the height of the
barrier was less than 1.0 kTr. The channel
currents obtained in the absence of any barrier and in the presence of a 1.0 kTr barrier were the same (23.7 ± 1.4 versus 22.1 ± 1.2 pA). With a further increase in the barrier
height, however, the current was systematically reduced, as shown in
Fig. 5. On the right-hand ordinate of
Fig. 5, we show the channel conductance at 150 mM concentration. It
decreased progressively from 110 pS with an increasing barrier height,
dropping to 52 ± 7 pS at 2.0 kTr. When the
barrier height was 2.5 kTr, this conductance of the channel became approximately the same as that obtained
experimentally from the ACh channel (Hamill and Sakmann, 1981) or from
the N-methyl-D-aspartate-activated channel
(Chung and Kuyucak, 1995).
|
Ionic concentrations in the channel
In the volume of our model channel, which is 2.16 × 1026 m3, there would be 720 water molecules.
At a concentration of 300 mM, a similar bulk volume would contain four
sodium and four chloride ions. Here we examine the number of sodium and
chloride ions inside the channel under various conditions. To compute
the average number and concentration of sodium and chloride ions inside
the channel, we divided the model catenary channel, whose length is 80 Å, into 16 5-Å sections, as shown in the inset of Fig.
6. The volumes of the slices from the
outermost layer to the smallest layer in the transmembrane segment were
2.94, 2.18, 1.88, 1.54, 1.16, 0.71, 0.25, and 0.15 × 1027 m3.
|
In the absence of a membrane potential and dipoles in the protein, ions were virtually excluded from entering the vestibule neck, owing to the repulsive dielectric force presented to them by the dielectric wall. In Fig. 6, A and B, the time averages of sodium and chloride concentrations in the channel are illustrated. The number of ions present in each layer per unit time was first tabulated (filled and open circles), and then the concentration in each layer was derived by taking into account its volume. The average ionic concentration in each of the two outermost layers was ~10% lower than that in the adjoining reservoir, which was 300 mM. The ionic concentration in successive layers declined progressively, dropping to less than 10% of the reservoir value in the neck region.
When a membrane potential of 100 mV was applied across the channel, such that the right-hand reservoir was made negative with respect to the left-hand side, there was a small but consistent increase in the concentration of sodium ions in the left vestibule and in the concentration of chloride ions in the right vestibule (Fig. 7). Other than this small asymmetry in the concentrations and a slight increase in the probability of ions being present in the constricted segment, the applied field had little effect on the average concentrations in the channel. A drastic change in the pattern of charge densities occurred when, instead of applying an electric potential across the channel, two rings of dipoles were placed on the channel wall. As shown in Fig. 8, there was a marked increase in the concentration of sodium ions in the constricted region of the channel. Sodium ions entering the channel occasionally would become detained in the potential well created by the dipoles on the wall. Some ions jumped from one well to the other and drifted across the other side, but because there was no potential gradient, the number of ions drifting in one direction (11.7 ± 1.3 pA) was about the same as that in the opposite direction (10.7 ± 1.3 pA).
|
|
To mimic the concentration gradient in the channel during its open state, we placed an energy barrier of 1.5 kTr at the constricted segment and applied a membrane potential of 100 mV. As shown in Fig. 9, the sodium concentration in all layers of the channel was approximately constant under these conditions and ions moved steadily from outside to inside, without being detained by the dipoles on the channel wall. In contrast, chloride ions were virtually excluded from entering the inside of the channel (Fig. 9 B). The conductance of the channel under this condition was 78 ± 4 pS at 150 mM (or 15.5 ± 0.7 pA at 300 mM), and no sodium ions traversed against the potential gradient.
|
Current-voltage relationships
The current-voltage relationships obtained from excised single channels are in general ohmic, although some show pronounced inward or outward rectification. The reversal potential observed in asymmetrical ionic solutions closely matches that predicted by the Nernst equation. Here we show how the presence of an energy barrier in the channel can distort the linear current-voltage relation.
The current increased linearly with the applied voltage, as shown in
Fig. 10 A, when there was no
additional potential barrier in the channel presented to ions for
penetrating the transmembrane segment. In this and all subsequent
current-voltage curves, we used the total strength of four dipoles in
each ring of 100 × 1030 Cm. The core conductance of
the channel, deduced from the regression line of the form,
I = V, fitted through the data points,
was 232 ± 4 pS (or 116 pS at 150 mM). When an energy barrier of
3.0 kTr was erected at the entrance of the
constricted segment, the current was attenuated, but not by a constant
proportion at all voltages. For example, the currents in the absence
and presence of this barrier at 200 mV were, respectively, 46.0 pA and
18.2 pA (39.6%). With the applied potential of 50 mV, however, the current was reduced from 14.0 pA to 2.2 pA (15.7%), thus indicating that the barrier of the same height became less of an impediment when
the driving force was large. With a 3.0 kTr
barrier, the current-voltage relationship became nonlinear, deviating
markedly from Ohm's law, as shown in Fig. 10 B. The solid
line fitted through the data points in Fig. 10 B was
calculated from the Pöschl-Teller function of the form
I = V/[1 + 
sech x],
where x = eV/VB and VB is the barrier height. The justification for
fitting the data with this function is given in the Discussion section.
The values of and used to generate the curves shown in solid
lines are given in the figure legend.
|
Fig. 11 A shows the current-voltage relationship obtained with asymmetrical ionic solutions in the two reservoirs. The ionic concentrations outside (left-hand reservoir) and inside (right-hand reservoir) were 480 mM and 120 mM, respectively. To make the solution electrically neutral, the same number of sodium and chloride ions was placed in each reservoir. With no energy barrier in the channel, the slope for the outward current (at the positive potential, filled circles) was steeper than that predicted by the Goldman equation. Then we added potassium ions such that the number of cations (sodium and potassium) in one reservoir was equal to that in the other reservoir. Thus the ionic concentrations outside were 480 mM NaCl and 120 mM KCl, whereas those inside were 120 mM NaCl and 480 mM KCl. Potassium ions were prevented from entering the constricted region of the channel by erecting two impermeable barriers at z = ±10 Å. Potassium ions colliding with these barriers were elastically scattered. The slope of the outward current (open circles) obtained under these conditions was less steep than that predicted by the Goldman equation. The solid line fitted through the data points was calculated with the Goldman equation of the form
|
(7) |
|
The slope of the outward current was affected more markedly by the presence of impermeable potassium ions when the concentration ratio was 2:1, as illustrated in Fig. 11 B. In these simulations, the ionic concentrations of NaCl outside and inside were, respectively, 400 mM and 200 mM. With no potassium ions, the magnitude of the inward current at a given applied potential is about half that of the outward current (filled circles). The currents flowing across the channel at various holding potentials were determined again with 200 mM KCl added to the extracellular solution and 400 mM KCl to the intracellular solution. Thus the total cation concentrations inside and outside in these simulations were 600 mM. As shown in Fig. 11 B, the outward current was appreciably attenuated, whereas the inward current remained unchanged (open circles). The solid line was again calculated from the Goldman equation, with the preexponential factor of 2 in the numerator of Eq. 7.
From a series of simulations such as those illustrated in Fig. 11, we conclude that the Nernst equation correctly predicts the reversal potential when the ionic concentrations in the two faces of the channel are different. The relative steepness of the slopes fitted through inward and outward currents changes appreciably when other cations that are impermeable to the channel are present. The shape of the current-voltage relationship is further distorted when an energy barrier is erected near the constricted segment of the channel.
Conductance-concentration curve
Experimentally, current across a biological ion channel increases
monotonically with an increasing ionic concentration initially and then
saturates with a further increase in concentration (Hille, 1992).
Saturation of channel currents occurs when there is a rate-limiting permeation process that is independent of ionic concentrations. For
example, an ion arriving near the constricted membrane segment will be
detained there for a period of time if, before traversing the narrow
pore, it needs to gain a sufficient kinetic energy to climb over an
energy barrier. The reason for the presence of such a barrier is
explained in the Methods.
In Fig. 12, the conductance of the
channel is plotted against the concentrations of sodium ions in the
reservoirs. The two reservoirs contained an equal number of
sodium-chloride pairs, and an applied membrane potential of 100 mV
provided the driving force for sodium ions to move inward. The presence
of two rings of dipoles, with their negative poles pointing to the
lumen, ensured that the channel was selectively permeable to sodium
ions. When the channel had no potential barrier, the ionic current
carried by sodium ions increased linearly with concentration, as shown in Fig. 12 A. Because the magnitude of the current in this
series of simulations was large, we used the total simulation period of
0.225 µs for each point shown in Fig. 12 A. The linear
conductance-concentration relation became distorted when a barrier was
placed 5 Å from each end of the cylindrical pore. Fig. 12 B
illustrates the conductance-concentration curve obtained from the
channel with a step potential barrier of 3.0 kTr. The ordinate of Fig. 12 B is
expanded, because the currents were greatly attenuated by the presence
of the barrier. At a low ionic concentration, the conductance was
nearly proportional to the ionic concentration. As the ionic
concentration was increased, however, the conductance increased less
with increasing concentrations. We fitted the points with the curve
calculated from the Michaelis-Menten equation of the form
![I=<FR><NU>I<SUB><UP>max</UP></SUB></NU><DE>1+K<SUB><UP>s</UP></SUB>/[c]</DE></FR>](/content/vol75/issue2/fulltext/793/img010.gif)
|
(8) |
|
| |
DISCUSSION |
|---|
|
Top
Abstract
Introduction
Methods
Results
Discussion
Concluding Remarks
References
|
|---|
As we have demonstrated here and elsewhere (Li et al., 1998), the
Brownian dynamics algorithm devised by van Gunsteren and Berendsen
(1982) is highly suitable for studying the motions of interacting ions
in a fluid. One of the major advantages of this algorithm is that the
range of time steps t is not restricted by the relaxation
time constant (1) as t 
1, which is the case in most other Brownian dynamics
algorithms. At room temperature, the relaxation time constant for
sodium ions is on the order of 10 fs; thus, to meet the condition
stipulated in this equation, a time step on the order of 1 fs must be
used with such algorithms. As the computational cost of covering a simulation period of a few microseconds would be prohibitively expensive, the use of such algorithms in studying the permeation of
ions across biological channels is severely limited. Our preliminary tests of the van Gunsteren-Berendsen algorithm in a simple periodic boundary and a toroidal dielectric boundary showed that many of the
important features of the motions of interacting ions in a fluid could
be captured even when a t as large as 100 fs is used. Among these important features are mean square displacement, the velocity distribution, and the conductance of bulk electrolyte solutions. Furthermore, 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.
The results of simulations are in agreement with experimental findings
in several respects. First, the channel conductance is close to that
determined experimentally for the ACh channel when the height of the
energy barrier is ~3.0 kTr or 1.23 × 1020 J. Second, the current-voltage relationship obtained
with symmetrical solutions is close to linear for moderate applied
potentials, as is the case with many biological channels. The
relationship obtained with asymmetrical solutions shows rectification,
with the reversal potential close to that predicted by the Nernst
equation. Third, the conductance-concentration curve has the same shape as those observed experimentally, although the saturation concentration is higher than expected. The model also makes a number of other predictions. To render the channel permeable to ions, several charge
groups must be placed in the protein wall to counteract the repulsive
dielectric force. An excessive number of charge residues, however,
allows sodium ions to flow against the potential gradient, so reducing
the conductance of the channel. Moreover, counterions do not screen
the image repulsion to any great extent, in either the vestibule or the
constricted region. Finally, current-voltage curves can be expected to
be nonlinear at large potentials, or for channels with large energy
barriers.
Suppression of currents by an energy barrier
A potential barrier presented to an ion, arising from the
dehydration and substitution process and the electrostatic interaction between ions and charge multipoles on the wall, will ensure that only
the ions with thermal energies E larger than the barrier height will be allowed to go through the neck. The probability 
of
an ion surmounting a barrier of height VB
follows from the Boltzmann distribution as
![𝒫<FENCE>E>V<SUB><UP>B</UP></SUB></FENCE>=<FR><NU>2</NU><DE><RAD><RCD>&pgr;</RCD></RAD></DE></FR> (kT)<SUP><UP>−</UP>3/2</SUP> <LIM><OP>∫</OP><LL><UP>V<SUB>B</SUB></UP></LL><UL>∞</UL></LIM> <UP>exp</UP>[<UP>−</UP>E/kT]<RAD><RCD>E</RCD></RAD> <UP>d</UP>E,](/content/vol75/issue2/fulltext/793/img011.gif)
|
(9) |
function,
|
(10) |
function,
Eq. 10 can be written as
|
(11) |
|
Current-voltage relationship
If the ionic concentrations on the two faces of the channel are the same, the current-voltage relationship obtained from patch-clamp recordings is, in general, ohmic. The current-voltage relationship deduced from our simulations becomes nonlinear whenever there is a potential barrier for the ions to surmount to cross the channel (Fig. 10). This deviation from Ohm's law is more pronounced when the potential barrier is large. Although the precise shape of the curve cannot be deduced a priori, it is easy to see how such a curvature in the current-voltage relationship would arise. The presence of a barrier is less of an impediment when the driving force is large. This intuitive observation suggests a modification of Ohm's law with the Pöschl-Teller function
). Then, the ionic strengths were changed
such that the outside contained 120 mM KCl (and 480 mM NaCl) and the
inside contained 480 mM KCl (and 120 mM NaCl), and the channel was made
impermeable to potassium ions. The outward current was markedly
depressed in the presence of potassium ions (
). (B)
The current-voltage relationship was obtained with ionic solutions of
400 mM outside and 200 mM inside (