A Combined Molecular Dynamics and Diffusion Model of Single Proton
Conduction through Gramicidin
Mark F. Schumaker,* Re´gis Pome`s,† and Benoît Roux‡
* Department of Pure and Applied Mathematics, Washington State University, Pullman, Washington 99164-3113 USA; † Theoretical
Biology and Biophysics Group, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 USA; and ‡ Groupe de Recherche en
Transport Membranaire, De´partements de Physique et de Chimie, Universite´ de Montre´al, Que´bec H3C 3J7, Canada
ABSTRACT We develop a model for proton conduction through gramicidin based on the molecular dynamics simulations
of Pome`s and Roux (Biophys. J. 72:A246, 1997). The transport of a single proton through the gramicidin pore is described
by a potential of mean force and diffusion coefficient obtained from the molecular dynamics. In addition, the model
incorporates the dynamics of a defect in the hydrogen bonding structure of pore waters without an excess proton. Proton
entrance and exit were not simulated by the molecular dynamics. The single proton conduction model includes a simple
representation of these processes that involves three free parameters. A reasonable value can be chosen for one of these,
and the other two can be optimized to yield a good fit to the proton conductance data of Eisenman et al. (1980, Ann. N.Y.
Acad. Sci. 339:8–20) for pH  1.7. A sensitivity analysis shows the significance of this fit.
GLOSSARY
Symbols that appear in two or more subsections are given.
When uppercase and lowercase symbols are given together,
lowercase denotes dimensionless quantities. Species s may
represent either H (proton) or d (defect). Roman R denotes
either side I or side II of the channel.
Latin symbols
CR Bulk concentration on side R (see Fig. 1 B)
Ds Effective standard diffusion coefficient for
species s (see Eq. 22)

s Diffusion coefficient of dipole moment reaction
coordinate (Eq. 17)

i Diffusion coefficient for independently tumbling
waters (Appendix A)
ˆ Laplace transform of diffusion coefficient (Eq.
16)
E Electric field in pore (Eq. 23)
e Base of natural logarithms
e0 Elementary electrical charge
G0 Single proton model conductance at
symmetrical equilibrium (see Eq. 21)
GV Chord conductance at potential V
Gmax Maximum apparent single proton conductance
(see above Eq. 21)
K Spring constant
KM Apparent Michaelis constant (see above Eq. 21)
k Rate of barrier crossing (Eq. 18 and Appendix B)
kB Boltzmann’s constant
I Current through channel (Eq. 21)
L Spatial length of channel (see Fig. 1 B)
ns Net number of waters in pore with dipole
moments oriented in z direction (Eqs. 2 and 7)
nmax Number of water molecules in pore (Eq. 5)
T Absolute temperature
ta Framework model free parameter (see
Optimization of Model Parameters in Methods)
t Time
t Mean first passage time
VI Applied transmembrane potential (see Fig. 1 B)
z Spatial coordinate coaxial with the pore (see
Fig. 1 B)
Greek symbols
 Friction coefficient (Appendix A)
 Framework model free parameter (see above Eq.
21)

s Dipole moment reaction coordinate for species s
(see Eqs. 3 and 8)
A
s Maximum extents of dipole moment intervals
(see Eqs. 6 and 10)

w Dipole moment of a water molecule (Eq. 1)

s Scaling factor between the orientation and
dipole moments for species s

s Potential of mean force for species s (see Fig. 2)

s Orientation moment for species s (see above
Eq. 1)
min
d Orientation moment coordinate of d potential
minimum
A
s Maximum extent of orientation moment interval
(see Fig. 2 and Eqs. 5 and 9)
Received for publication 14 February 2000 and in final form 1 September
2000.
Dr. Pome`s’ present address is Department of Structural Biology and
Biochemistry, Hospital for Sick Children, Toronto, ON M5G 1X8, Canada.
Dr. Roux’s present address is Department of Biochemistry and Structural
Biology, Weill Medical College of Cornell University, 1300 York Ave.,
New York, NY 10021.
Address reprint requests to Dr. Mark F. Schumaker, Department of Pure
and Applied Mathematics, Washington State University, Pullman, WA
99164. Tel.: 509-335-7273; Fax: 509-335-1188; E-mail: schumaker@
wsu.edu.
© 2000 by the Biophysical Society
0006-3495/00/12/2840/18 $2.00
2840 Biophysical Journal Volume 79 December 2000 2840–2857

B
d Effective electrical coordinates of boundary
regions (see Fig. 2 B)
C
d Maximum extent of interior of defect interval
(see Fig. 2 B)

w Orientation moment of a water molecule (Eq. 1)
 Electrical potential energy of pore contents (Eq.
23)
INTRODUCTION
Proton conduction through gramicidin has gained recent
attention because 1) Single-channel currents are large and
can be measured over a very wide range of concentrations.
The observed conductance has a very prominent shoulder as
a function of [H] between pH 1 and pH 2 (Eisenman et al.,
1980). This observation has been confirmed by Akeson and
Deamer (1991), and a similar shoulder has recently been
observed in the RR dioxolane gramicidin analog by Cu-
FIGURE 1 (A) Orientation of water molecules. Water molecules tend to
be oriented so that each shares one hydrogen bond with the channel and
two with neighboring waters. The axial component of the water dipole
moments may reverse about a water that shares two bonds with the
channel. This 2D depiction of hydrogen bonds is stylized; no true asym-
metry about the pore axis is implied. (B) Applied potential. The single-
proton model assumes that the electrical potential, V, is a linear function of
the coordinate parallel to the pore axis, z. Side I refers to the solution on the
left, and side II refers to the solution on the right. VI is the applied electrical
potential on side I. VII  0 by convention. (C) Hypothetical single-proton
conduction mechanism. The top row of cartoons shows a pore with an
excess proton. The proton may enter the channel on side I at the upper left.
In this cartoon, the orientation of the pore waters favors entrance on side
I. As the proton passes through the pore, the hydrogen bond structure of
surrounding waters is altered so that water oxygens remain oriented toward
the center of excess charge. The bottom row shows a pore in the absence
of an excess proton. Water dipole moments flip as a defect in the hydrogen
bonding structure diffuses through the column of waters. (D) State diagram
of the single-proton model. The top segment corresponds to the proton-
occupied pore and is parameterized by the coordinate H. The bottom
segment corresponds to a defect traversing the central barrier in Fig. 2 B
and is parameterized by d  [C
d , C
d ]. The endpoints of the bottom
segment are the boundary states bI and bII. They correspond to a defect in
the boundary regions C
d
 
d
  A
d .
FIGURE 2 Proton and defect potentials of mean force. All energies are
given in units of kBT for T  298 K. (A) The ordinate is the proton
potential,H, in kBT, and the abscissa is the orientation moment of the pore
contents, H. Dipole moments are obtained by scaling H  H H. (B)
The ordinate is the defect potential, d, in kBT, and the abscissa is the
orientation moment, d. Protons may enter the channel when defects are in
the boundary regions I and II. Shown in the figure are wide boundary
regions defined by A
d
 8.1, B
d
 6.8, and C
d
 5.7 e0Å. The text also
considers narrow boundary regions with A
d
 8.1, B
d
 7.44, and C
d

6.9 e0Å. Boundary regions surround the interior of the defect interval,
defined by C
d

d
C
d . Dipole moments are obtained by scaling d 

d

d.
Single H Conduction through Gramicidin 2841
Biophysical Journal 79(6) 2840–2857

kierman (2000). 2) Experimental evidence suggests that the
proton permeation mechanism is qualitatively different
from that for other cations (Busath et al., 1998; Phillips et
al., 1999). Molecular dynamics simulations suggest that
water reorientation may be involved in the rate-limiting
steps (Pome`s and Roux, 1997, and manuscript in prepara-
tion). 3) Proton conduction can be used to probe permeation
by other cations (Heinemann and Sigworth, 1989). 4) Pro-
ton conduction through gramicidin provides an interesting
comparison with other proton channels (DeCoursey and
Cherney, 2000). Among the latter are channels that form
parts of molecular machines that convert proton gradients
into chemical energy by synthesizing ATP (for example,
Boyer, 1997) or convert proton gradients into mechanical
energy by turning flagellar rotors (for example, Elston and
Oster, 1997).
This paper develops a model for proton conduction
through gramicidin based on the molecular dynamics sim-
ulations of Pome`s and Roux (1997 and manuscript in prep-
aration). These simulations incorporate the detailed struc-
ture of the channel (Arseniev et al., 1985). Techniques
developed to describe the dynamics of this relatively simple
molecule may soon be applied to channels of broader bio-
logical importance. Nuclear magnetic resonance and x-ray
crystallography have recently provided high-resolution
structures for several ion channels (Cowan et al., 1992;
Ketchem et al., 1993; Doyle et al., 1998; Chang et al.,
1998).
A large gap separates the time scales of molecular mo-
tions described by molecular dynamics and the time scales
required to form meaningful averages of ion permeation
events comparable to experimentally observed conduc-
tances. In a companion paper (Schumaker et al., 2000b) we
develop a framework model for ion permeation that bridges
the gap. This probabilistic model is designed to combine the
results of molecular dynamics with simple assumptions
about the proton entrance and exit processes that were not
addressed by the simulations. The framework model is
based on two coupled Nernst-Planck equations, one describ-
ing diffusion of a single proton and the other describing
diffusion of a defect in the hydrogen bond structure of pore
waters. The interaction between proton permeation and the
dynamics of the pore waters can be understood in terms of
a state diagram analogous to those encountered in rate
theory.
The present paper is chiefly concerned with incorporating
the results from molecular dynamics into the framework
model and making a detailed comparison with experiment.
We discuss the proton and water reorientation potentials of
mean force and deconvolute their velocity autocorrelation
functions to obtain numerical values for diffusion coeffi-
cients. We also verify that a mathematical approximation
made by the framework model, the lumped state approxi-
mation, is accurate when applied to the water reorientation
potential of mean force. After incorporating the information
from molecular dynamics, the model has three free param-
eters. These are concerned with proton entrance and exit—
processes not addressed by the simulations.
Values for these parameters are obtained from the con-
ductance and conductance ratio data of Eisenman et al.
(1980), who present an extensive set of measurements for
pH  1.7, where we believe single proton conductance
dominates. A reasonable value is chosen for one parameter,
and the other two are optimized to fit the data. A good fit is
obtained, and the resulting model is also compared with two
I–V curves measured by Decker and Levitt (1988). We also
discuss our model in light of the recent experimental results
of Phillips et al. (1999) on proton conduction in gramicidin
analogs.
METHODS
Dipole moment of the gramidicin pore contents
The pore of gramicidin is roughly cylindrical, and conse-
quently an experimentally applied transmembrane potential
gives rise to an approximately constant electric field within
the pore (Jordan, 1982; Jordan et al., 1989; Roux, 1999). A
sketch of the pore geometry is given in Fig. 1 A. When the
electric field is constant, the component of the pore potential
energy due to the applied field depends only on the net
charge and the z component of the dipole moment of the
pore contents (Schumaker et al., 2000b). In part for this
reason, we use the z component of the dipole moment of the
pore contents as the reaction coordinate for the single-
proton model.
The reaction coordinate of the molecular dynamics sim-
ulation of Pome`s and Roux 1997; manuscript in prepara-
tion) depends on the orientation moment of the polarizable
and dissociable PM6 waters (Stillinger and David, 1978;
Weber and Stillinger, 1982) occupying the pore. Fig. 2, A
and B, shows the potential of mean force of the pore
contents with and without an excess proton, respectively, as
a function of the z component of the orientation moment.
This section describes what the orientation moment is and
how it is rescaled to obtain an estimated dipole moment of
the water column. We refer to a channel pore without an
excess proton as an empty pore.
The orientation moment of a single water molecule is the
formal dipole moment obtained by assigning oxygen a
formal charge of 2e0 and hydrogen a formal charge of
1e0. The OH separation of the model water is 1 Å, and
the H-O-H angle is 110°. Thus the orientation moment of
the PM6 water molecule is 2 cos 55°
1.15e0 Å. Based
on the work of Duca and Jordan (1997; see Results), we use
the constant value of 2.4 Debye  0.5e0 Å for the water
dipole moment in both empty and occupied pores. We
therefore introduce a scaling factor of d  0.5/1.15 be-
tween the orientation moment and the estimated water di-
pole moment.
2842 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

A water typically forms one hydrogen bond with a pore
backbone carbonyl and additional hydrogen bonds with the
neighboring waters on either side. This configuration is
suggested in Fig. 1, A and C, although in a highly stylized
representation. The water dipole moments are tilted from
the z axis by an average angle of 45°. We use the value

w
 (1.15e0 Å) cos 45
0.81 e0Å for the z component of
the water orientation moment. Then the z component of the
water dipole moment, w, is given by

w


d

w . (1)
The reaction coordinate d of an empty pore is the z
component of the dipole moment of the pore waters. d is
the z component of the orientation moment of the pore
waters. These coordinates parametrize the progress of a
defect in the hydrogen bond structure of the pore waters
(Pome`s and Roux, manuscript in preparation). Let nd be the
net number of pore waters with dipole moments oriented in
the z direction: the number of waters with z orientation
minus the number with z orientation. d and d can be
expressed in terms of nd:

d

ndw , (2)

d

ndw . (3)
Then the relationship between the dipole and orientation
moments of an empty channel is simply

d


d

d . (4)
Let nmax be the number of pore waters. To define the
single file based on simulations, one can either count only
those water molecules that can make up to two hydrogen
bonds with other water molecules, or count only waters that
are confined near the central axis of the channel pore.
nmax  10 is a consistent choice according to both criteria
(Pome`s and Roux, 1996). We must have nmax n
d

nmax. Let A
d and A
d be the maximum magnitudes of d and

d, respectively, corresponding to configurations with all
waters aligned. We then have A
d

d
A
d and A
d


d
A
d , where
A
d

nmaxw , (5)
A
d

nmaxw , (6)
and A
d
 
d
A
d . With the value of w obtained above, we
have A
d

8.1 e0Å. In Fig. 2 B this value corresponds to the
inflection points exterior to the potential minima on either
side of the central barrier. Although the range of d values
formally extends to 9.1 e0Å, the most extreme values
correspond to strained configurations that are probably of
little physical importance.
When the pore is occupied by an excess proton, the dipole
moment includes a contribution due to the excess charge
calculated about the center of the pore, z  0:

H

e0z nHw , (7)

H

e0z nHw , (8)
where nH is the net number of pore waters with dipole
moments oriented in the z direction in the occupied chan-
nel. We assume nmax n
H
nmax. If the pore occupies
the interval L/2 z L/2, it follows that A
H

H

A
H and A
H

H
A
H, where
A
H

e0L/2 nmaxw , (9)
A
H

e0L/2 nmaxw . (10)
In the occupied channel, the idealized state of maximum
dipole moment corresponds to the excess charge at the
extreme coordinate z  L/2, with nmax pore water dipole
moments aligned in the z direction (see the upper right
cartoon of Fig. 1 C). In this way, the partial negative charges
carried by water oxygens are closest to the center of positive
charge.
Combining Eqs. 5 and 9 and Eqs. 6 and 10, we have
A
d
 A
H

e0L/2
A
d
A
H . (11)
From the molecular dynamics, A
H
 3.35 e0Å; then using
A
d
 8.1 e0Å, we obtain L  22.9 Å, 2 Å shorter than the
physical length of the pore. This discrepancy may be ex-
plained, at least in part, by the procedure used to simulate
the channel occupied by an excess proton. In these simula-
tions, the pore is occupied by 10 PM6 water molecules, and
there are also several TIP3P waters clustered about the
channel entrances at either end. The TIP3P waters cannot
form hydrogen bonds. The closest that the excess proton can
approach the channel entrances corresponds to hydronium-
like (H3O
) configurations on the outermost PM6 waters. If
the TIP3P waters were able to host hydrogen bonds, there
would be O2H5
 ions between the PM6 and TIP3P waters.
Therefore, the value of A
H obtained from the simulations
may be slightly low, and, as a result, the effective length of
the simulated channel may be slightly shorter than the pore
of gramicidin.
Finally, we will obtain a scaling factor relating the mo-
lecular dynamics reaction coordinate H and the pore dipole
moment H. We make the simplifying assumption that
when the excess proton is located at coordinate z, all water
dipole moments will be oriented away from it, with water
oxygens aligned toward the excess positive charge. Then nH
will be a linear function of z. Put
nH
 2/L
nmaxz , (12)
so that when z  L/2 (the proton is at the entrance on side
II), nH  nmax. Substituting Eq. 12 into Eqs. 7 and 8, and
then dividing, we have

H


H

H , (13)
Single H Conduction through Gramicidin 2843
Biophysical Journal 79(6) 2840–2857

where

H

e0L 2nmaxw
/e0L 2nmaxw
. (14)
Diffusion coefficients
Diffusion coefficients were calculated following the method
of Crouzy et al. (1994) and Woolf and Roux (1994), based
on the Generalized Langevin equation (Berne et al., 1988;
Straub and Berne, 1988; Berne and Pecora, 1976):
m¨ 

0
t
Mt

˙

d
 K
Ft
, (15)
where  is the reaction coordinate (the z component of the
dipole moment of the pore content), m is a reduced mass, K
is a spring constant, M is the memory kernel (a generalized
friction coefficient), and a dot denotes the derivative with
respect to time. F(t) represents the influence of degrees of
freedom orthogonal to the reaction coordinate. Assuming
that these other degrees of freedom relax on much shorter
time scales than the reaction coordinate, F(t) is modeled as
a random process.
If Eq. 15 is multiplied by ˙(0) and an ensemble average
is taken, an equation is obtained for the velocity autocorre-
lation function of . The result is then Laplace transformed,
and the Einstein relation is applied to find a formula for the
Laplace transform of the diffusion coefficient (Crouzy et al.,
1994; Woolf and Roux, 1994):
ˆs

 cˆs
2˙2
cˆs
s2 ˙2/s 2˙2
. (16)
In this expression, cˆ(s) is the Laplace transform of the
velocity autocorrelation function c(t)  ˙(t)˙(0), and
. . . denotes ensemble average. The diffusion coefficient is
obtained as the limit

lim
s3 0
ˆs
. (17)
For simplicity, we have assumed that the proton and defect
diffusion coefficients are constant, independent of the value
of their reaction coordinates. In both cases, c(t) was com-
puted from molecular dynamics simulations at 0.5-fs inter-
vals for times in the range 0  t  0.5 ps. Fig. 3 A shows
the normalized proton velocity autocorrelation function

H(t)H(0) for the proton reaction coordinate held near

H
 2.99 e0Å by an umbrella potential. Values of the
proton orientation moment are scaled to the dipole moment
by H  HH, as described in the previous section. Fig. 4
A shows ˆH as a function of s. Parabolic extrapolation of
ˆ
H from the interval 10  s  40 gives an estimate H 
3.52  (H)2 (e0 Å)
2 ps1. The smooth behavior of ˆH at
high values of s becomes corrupted for s  10 correspond-
ing to times t  0.1 ps. This may reflect other processes
influencing the transport of protons on longer time scales,
for example, interactions between proton motion and rear-
rangement of hydrogen bonds.
Similarly, the defect diffusion coefficient was estimated
from the defect reaction coordinate held near d 
0.104 e0Å by an umbrella potential. Values of the defect
dipole moment are given by d  dd. Fig. 3 B shows the
normalized defect velocity autocorrelation function. Fig. 4
B shows ˆd as a function of s. Parabolic extrapolation of ˆd
from the interval 4  s  10 gives an estimate d  1.08 
(d)2 (e0 Å)
2 ps1. Values of ˆd become corrupted for s 
3 corresponding to times t  0.33 ps; this reflects the finite
interval of times for which c(t) was computed.
The lumped state approximation
The framework model lumps the boundary regions of the
defect potential of mean force, shown in Fig. 2 B, into
discrete boundary states. These are represented by the filled
endpoints of the bottom segment of the state diagram shown
by Fig. 1 D. The weights of these boundary states are
FIGURE 3 Velocity autocorrelation functions. (A) The normalized pro-
ton velocity autocorrelation function H(t)H(0) for the reaction coordi-
nate held near H  2.99 e0Å by an umbrella potential. The unnormalized
amplitude is H (t)2  7346.6 (e0 Å)
2 ps1. (B) The normalized defect
velocity autocorrelation function for the reaction coordinate held near d 
0.104 e0Å by an umbrella potential. The unnormalized amplitude is

d(t)2  5270.8 (e0 Å)
2 ps1.
2844 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

carefully adjusted to give the integral over the Boltzmann
distribution under the conditions of a symmetrical equilib-
rium (Schumaker et al., 2000b). This lumped state approx-
imation (LSA) makes it possible to solve the framework
model analytically. However, the boundary states act elec-
trically like points with electrical coordinates B
d , where
B
d is a free parameter of the framework model. The value of
B
d is chosen to give the best possible representation of the
mean first passage time (MFPT) over the defect potential
barrier as a function of the applied potential VI.
To determine the best possible representation of the
MFPT, consider a diffusion process that does not depend on
the lumped state approximation. This latter process is de-
fined on the interval [A
d , A
d ] of the potential profile
shown in Fig. 2 B. The process starts at d min
d
6.5
e0Å, a minimum of 
d, with a reflecting barrier at d 
A
d . The mean first passage time to d  min
d is then
calculated as described in Appendix B. Because the MFPT
found in this way is the standard for comparison, we call it
the exact MFPT. An analytical formula can also be obtained
for LSA MFPTs, as will be described by Mapes and Schu-
maker (manuscript submitted for publication). These mean
times correspond to a process that starts at the lumped state
at d  C
d . The exact mean first passage time to the
lumped state at C
d is then given. Values of the electrical
coordinate for the lumped states, B
d , are chosen so that LSA
MFPTs approximate the exact values well over the range of
applied potentials, 200 mV VI 200 mV.
Fig. 5 A compares the LSA MFPTs over the central
barrier shown in Fig. 2 B with the exact MFPTs. The
logarithm of the MFPTs is shown as a function of the
applied potential VI. When VI  0, the potential energy of
the well centered at d  min
d is increased relative to the
energy of the well at min
d , reducing the MFPT. To test the
sensitivity of our results to the value of C
d , MFPTs for the
LSA are calculated for both wide and narrow boundary
regions. The wide boundary regions are shown in Fig. 2 B
FIGURE 4 Estimates of diffusion coefficients. In each figure, the nu-
merical estimate of the Laplace transform of the diffusion coefficient, ˆ, is
shown as a function of the transform parameter s. (A) Values of ˆH(s).
Parabolic extrapolation to the ordinate gives the estimated proton diffusion
coefficient of H  3.52  (H)2 (e0 Å)
2 ps1. The Laplace transform is
shown for H  1. (B) Values of ˆd(s). The estimated defect diffusion
coefficient is d  1.08  (d)2 (e0 Å)
2 ps1. The transform is shown for

d
 1.
FIGURE 5 Logarithm of the mean first passage time over the defect
potential central barrier plotted as a function of the applied electrical
potential. Squares and circles give estimates of the mean first passage time
that apply to the potential defined in the full interval betweenA
d . Squares
give the exact mean first passage between the potential minima, with
reflecting boundary conditions at A
d . Circles give the Kramers approx-
imation to the mean first passage time. Triangles give mean first passage
times, using boundary states with C
d
 5.7 e0Å. Diamonds give times
using boundary states with C
d
 6.9 e0Å. (A) Defect potential given by
molecular dynamics. (B) Defect potential rescaled so that the height of the
central barrier is 4kBT lower than that given by molecular dynamics.
Single H Conduction through Gramicidin 2845
Biophysical Journal 79(6) 2840–2857

and have C
d
 5.7, B
d
 6.8, and A
d
 8.1 e0Å. The narrow
boundary regions have half the width of the wide regions,
with C
d
 6.9, B
d
 7.44 e0Å, and the same value of A
d .
Also shown in Fig. 5 A are MFPTs estimated by using
Kramers’ escape theory (e.g., Risken, 1989). The integrals
in the Kramers formula are calculated numerically so that a
robust estimate of the MFPT is obtained, even though the
crest of the central barrier in Fig. 2 B has a complex shape;
this is described in Appendix B. As shown by Fig. 5 A,
MFPTs obtained by all four methods are in very good
agreement over the full range of applied potentials consid-
ered.
The sensitivity analysis presented in the Results depends
on the LSA remaining accurate for a range of barrier
heights. In one case, the amplitude of the potential profile
shown in Fig. 2 B is reduced by 4kBT, leaving a barrier that
is little more than 2kBT high in the absence of an applied
potential. Fig. 5 B shows how the LSA MFPTs then com-
pare with the other MFPT values. Both sets of LSA MFPTs
remain in good agreement with the exact values. However,
estimates of the MFPT obtained from the Kramers formula
diverge from the other values as VI increases. This is not
surprising, because Kramers’ theory is an asymptotic ap-
proximation that is accurate when potential barriers are
higher than 4kBT (Cooper et al., 1985, 1988). However,
when VI  200 mV, the height of the potential barrier as
seen by a diffuser beginning on side I is less than 1kBT.
Fig. 5 shows that the logarithm of the mean first passage
time is approximately a linear function of VI. An interesting
comparison can be made between linear least-squares fits to
the points in the figure and the expected expression in the
limit of high potential barriers. According to Kramers’
theory, the rate k of barrier crossing in the latter limit has the
form
k
Aexp Wd/kBT , (18)
where Wd is the height of the total defect potential barrier
and A is a constant. The total potential, Wd, is the sum of the
potential of mean force, d, and the electrical potential
energy, d (Schumaker et al., 2000a). The mean first pas-
sage time corresponding to rate k is 1/k. The height, Wd,
depends on the applied potential according to
Wd
W0
d
 fe0VI , (19)
where W0
d is the height in the absence of an applied
electrical potential and f is the fraction of the electrical
potential drop between the maximum and minima of d. To
be precise, put f  dmin
d /(e0L), where min
d
 6.5 e0Å is
the coordinate of the defect potential minimum and e0L 
22.9 e0Å is the total potential drop associated with the
passage of a proton across the membrane (see Eq. 11). The
predicted slope of log10 (mean first passage time) as a
function of VI is then m, where
m
log10 e
fe0/1000kBT
, (20)
and where e is the base of the natural logarithms, e0 is the
elementary electrical charge, and the factor of 1000 is due to
measuring VI in mV.
Using the value of d described above (in Dipole Mo-
ment of the Gramicidin Pore Contents), we obtain, from Eq.
20, m  0.00211. In comparison, linear least-squares fits
through each of the four sets of values in Fig. 5 A give
slopes in the range 0.00214 m 0.00216, close to the
limiting value. However, linear least-squares fits through
the exact and LSA sets in Fig. 5 B give 0.00189 m
0.00191. For the W0
d

2kBT potential barrier used to
generate the latter figure, the LSA gives a softened depen-
dence of mean first passage time on the applied field, which
is similar to the exact result. In contrast, the slope obtained
from the Kramers’ estimates in Fig. 5 B gives m  0.00209.
Optimization of model parameters
The framework model of proton conduction (Schumaker et
al., 2000b) has three free parameters: , ta, and C
d . 
depends on the absolute free energy difference between the
proton and defect potentials of mean force, and ta controls
the rate of proton exit from the channel; together  and ta
control the rate of proton entrance and exit. C
d determines
the width of the boundary regions of the defect potential
profile. These are defined as regions within which there may
be a defect in the hydrogen bond structure of waters when
a proton enters or leaves the pore (see Fig. 2 B).
Initial estimates of  and ta are based on subjective
estimates of the Michaelis constant, KM, and the maximum
conductance, Gmax, associated with entry of the first proton
into the gramicidin channel. KM and Gmax are assumed to be
associated with the lower leg of rising conductance in the
data of Eisenman et al. (1980), which is measured at sym-
metrical pH and a transmembrane potential of 50 mV (see
Fig. 7 A). The estimated values are KM  0.005 M and
Gmax  24 pS.
We use these estimates to determine values of X and Y in
the formula for G0, the conductance at symmetrical concen-
trations and zero transmembrane potential (Schumaker et
al., 2000b):
G0

I
VI

0


e0
2
kBT
C
X CY C2Z
. (21)
In this formula, I is the proton current and VI is the trans-
membrane potential. The subscript 0 refers to the derivative
evaluated at the equilibrium state with symmetrical concen-
trations CI  CII  C and VI  0. The coefficients X, Y, and
Z depend on the three free parameters in the following way:
X  X(ta, ), Y  Y (ta, C
d ), and Z  Z (, C
d ).
Assuming a fixed value of C
d , the estimates of X and Y
determine initial values of ta and . This is used as a starting
point for local least-squares minimization of the model
against the combined data in Fig. 7, A and B, for pH  1.7.
2846 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

ta and  are varied to minimize the sum of unweighted
differences between the log10G50 data (Fig. 7 A) and the
model plus twice the sum of the unweighted differences
between the conductance ratio data (Fig. 7 B) and the
model. We consider two different values of the third pa-
rameter controlling the width of the boundary regions: C
d

5.7 e0Å, giving wide boundary regions, and C
d
 6.9 e0Å,
giving narrow boundary regions. These are used to generate
the fits shown in Fig. 7, A and B.
RESULTS
Molecular dynamics simulations
Pome`s and Roux (manuscript in preparation) have simu-
lated the dynamics of proton conduction through gramici-
din. The orientation of the channel pore in the membrane is
sketched in Fig. 1 A. The packing of the pore contents in
several configurations is represented in simplified form by
Fig. 1 C. In one set of simulations the pore includes an
excess proton; these are represented by the top row of the
figure. The water oxygens tend to be oriented toward the
center of excess charge. The resulting water dipole moments
partially offset the contribution of the excess proton to the
total pore dipole moment.
In a second set of simulations, there is no excess proton.
These are represented by the bottom row of Fig. 1 C. It is
energetically favorable for interior pore waters to share one
hydrogen bond with the channel and two others with the
adjacent water molecules. In this way they tend to have their
dipole moments aligned with similar z components. How-
ever, pore waters near the channel entrances tend to be less
ordered. As a result, there is often a localized disorder in the
packing structure of the pore waters near the entrance,
centered on a water molecule that shares two hydrogen
bonds with the channel. This disordered region may diffuse
from one end of the pore to the other, reversing the dipole
moments of the pore waters.
Potentials of mean force
Fig. 2 shows the potentials of mean force plotted as a
function of the orientation moment H or d, the reaction
coordinate of the simulations. In Methods, we discuss how
these orientation moments can be scaled to obtain estimated
components of the dipole moment of the pore contents
parallel to the pore axis: d  dd and H  HH. Dipole
moments are calculated with respect to the center of the
pore, at z  0. According to the notation introduced in Fig.
2, axial components of the proton dipole moments vary over
the range A
H

H
A
H, and defect components vary
over A
d

d
A
d . The total interval of dipole moment
spanned by the two potentials is 22.9 e0Å, corresponding to
an elementary charge passing through the length of a pore
22.9 Å long. This is somewhat shorter than the physical
length of the gramicidin pore, which is 25 Å. A possible
explanation for this discrepancy is offered in Methods.
The scaling factor d is calculated by assuming that the
average dipole moment of a water molecule in the pore is
2.4 Debye. By comparison, Duca and Jordan (1997) have
estimated effective average dipole moments of water mol-
ecules in a simulated polyglycine analog of gramicidin. In
their simulations with all groups polarizable, the average
water dipole moment in an empty pore was found to be
2.3 Debye. In a pore occupied by a single Na or Cs ion,
this moment depends on the proximity of the water to the
ion and varies between 2.3 Debye far from the ion (e.g.,
fifth nearest neighbors) to 3.0 Debye next to the ion.
The value of A
d used corresponds to the z component of
the dipole moment of a column of 10 aligned water mole-
cules, each tilted at an angle of 45° from the pore axis,
consistent with the conformation of waters in the pore.
FIGURE 6 Single-proton model trajectories. Trajectories for A and B are
obtained by simulating a random walk approximating the single-proton
model for C
d
 5.7 and A
d
 9.1e0 Å, CI  CII  0.01 M, and VI  0.
(A) A proton trajectory. (B) A defect trajectory. Occupation of the bound-
ary states is indicated at C
d
 5.7 e0Å. (C) Defect trajectories diffusing
in the full interval of reaction coordinate defined by the molecular dynam-
ics simulations.
Single H Conduction through Gramicidin 2847
Biophysical Journal 79(6) 2840–2857

Calculation of the scaling factor H also assumes that the
center of excess charge in the occupied pore separates two
subcolumns of water molecules, both with dipole moments
aligned so that oxygens are brought close to the excess
proton. Details of these calculations are presented in Methods.
The proton potential is a shallow well. As a consequence,
charge density corresponding to an excess proton concen-
trates near the center of the pore. The defect potential has
two minima separated by an 6kBT barrier. Thus the defect
concentrates in regions near the channel entrances, on either
side of the central barrier. The fine structure of the defect
potential is due to both the periodicity of the gramicidin
helix and the energetics of breaking hydrogen bonds be-
tween the single-file water molecules. The absolute free
energy difference between the proton and defect potentials
is not known.
Ideally, the potential of mean force should include elec-
trostatic interactions with the membrane and bulk electro-
lyte (Roux, 1999). Instead, the potentials H and d depend
only on the channel, pore waters, and several additional
waters outside the channel entrance. Below, we will present
a sensitivity analysis showing how our conclusions change
as we perturb the defect barrier height and diffusion
coefficient.
Diffusion coefficients
Diffusion coefficients of protons and defects simulated by
molecular dynamics were calculated by the Laplace Trans-
form method (Woolf and Roux, 1994; Crouzy et al., 1994),
as described in the Methods. We obtained H  3.52 
(H)2 (e0Å)
2 ps1 and d  1.08  (d)2 (e0Å)
2 ps1

0.20 (e0Å)
2 ps1. The defect value can be compared with a
simple estimate of the dipole diffusion coefficient i for a
model of 10 independently tumbling water molecules. In
Appendix A we obtain i  0.36 (e0 Å)
2 ps1, comparable
to the values of d given above.
The value of H can be associated with a standard
diffusion coefficient in the following way. Consider a pro-
ton starting at one channel entrance and diffusing to the
other entrance. In the absence of a potential and assuming a
and diffusion coefficients, and the wide boundary regions shown in Fig. 2
B; optimized values are ta  23.5 ns and   3.60. Long-dashed curves
show fits using a defect potential of mean force scaled to an amplitude
2kBT higher than that given by molecular dynamics; the proton potential of
mean force and all diffusion coefficients are those given by molecular
dynamics. Optimized values are ta  16.5 ns and   3.53. Short-dashed
curves show fits using molecular dynamics estimates of potentials and
diffusion coefficients and narrow boundary regions. Optimized values are
ta  21.6 ns and   2.23. (A) Fits to conductance measured at 50 mV
applied potential. (B) Fits to conductance ratios. The experimental data are
conductances measured at 100 mV divided by estimated conductances at 0
mV. (C) Comparison between framework model fits to data of Eisenman
et al. and the experimental I–V curve of Decker and Levitt (1988) at pH
3.75. (D) Comparison with Decker and Levitt I–V curve at pH 2.75.
FIGURE 7 Comparisons between the framework model and experimen-
tal data. Circles are the data of Eisenman et al. (1980). Curves show
framework model fits to data. Fits are made optimizing the values of ta and
, parameters that control the rate of proton entry and exit. Solid curves
show fits using molecular dynamics estimates of potentials of mean force
2848 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

reflecting boundary at the starting point, the mean first
passage time to reach the goal is t  (2A
H)2/(2H). The
corresponding mean first passage time to diffuse the phys-
ical length, L, of the pore with a standard diffusion coeffi-
cient DH is t  L2/(2DH). Equating these two expressions,
we obtain
DH
HL2/2A
H

2 . (22)
Since A
H
 
H
A
H, factors of H cancel in this formula.
Using the value of A
H
 3.35 e0Å from the molecular
dynamics and channel length L  22.9 Å, we obtain DH

41 Å2 ps1. This is 40 times the diffusion coefficient of
protons in water, 9.3  105 cm2 s1  0.93 Å2 ps1
(Hille, 1992), but only about twice as great as the diffusion
coefficient of protons in ice (Eisenberg and Kauzmann,
1969).
The Laplace Transform method of estimating the diffu-
sion coefficient requires extrapolation from a smooth curve,
which is determined by an analysis over short time scales. In
the case of Fig. 4 A, extrapolations are made from values of
s corresponding to time scales of 0.1 ps. Mechanisms
operating on longer time scales, such as movement of de-
fects in the water chain, may decrease the effective diffusion
coefficient.
Single-proton conduction model
Hypothetical conduction mechanism
We will construct a kinetic model based on the assumption
that protons can only enter the channel under the energeti-
cally favorable circumstance that the net dipole moment of
the pore contents is pointing away from the proton. The
results obtained from the simulations then suggest a mech-
anism of proton conduction through gramicidin.
The geometry of the pore between aqueous solutions on
sides I and II is indicated by Fig. 1 A. Sides I and II are also
indicated for the upper left cartoon in Fig. 1 C. As suggested
by Fig. 1 C, proton entry into the channel from side I would
be facilitated by waters lined up to present a partial negative
charge at the channel entrance. As the proton passes through
the channel, waters remain oriented to the center of excess
charge. When the proton escapes the dipoles remain nearly
aligned. A defect must then propagate through the pore to
restore the original orientation of waters in the channel. An
elementary charge is transferred from side I to side II when
the channel cycles once around the diagram clockwise.
The lower middle cartoon in Fig. 1 C depicts the defect
separating two oppositely aligned columns of water. The
dipole moments of both columns point toward the defect,
giving it a partial positive charge. This is the nature of the
defect most commonly encountered in the simulations. It is
called entrance-initiated by Phillips et al. (1999), because it
originates near the channel entrance opposite the exiting
proton. Alternatively, the dipole moments of the two partial
columns may point away from the defect, giving it a partial
negative charge (an exit-initiated defect). The mathematical
formulation of the single-proton model does not depend on
this choice. Our use of the word “defect” here is slightly
different from that encountered in reference to Bjerrum L or
D defects (for example, Eisenberg and Kauzmann, 1969). In
this latter case, a defect refers to an interruption of the
hydrogen bond chain. The terminology of Phillips et al. has
the advantage of avoiding confusion, because an entrance-
initiated defect is likely to propagate through the formation
of transient Bjerrum L defects (Pome`s, 1999), which are
sometimes referred to as negative defects.
Applied transmembrane potential
Fig. 1 B shows the electrical potential, V(z), due to a voltage
difference, VI, between the two sides of the membrane.
Suppose VI  0. In the top row of Fig. 1 C, the field pushes
the proton toward side II. However, the dipole moment of
the pore contents changes sign once the proton leaves the
channel, and so the electric field then tends to reverse the
orientation of the waters, now pushing the defect from side
I to side II. A positive potential on side I favors clockwise
cycling around the diagram.
For the relatively uniform cylindrical geometry of the
gramicidin pore, the applied electrical potential is well
approximated as decreasing linearly across the length of the
pore (Jordan, 1982; Jordan et al., 1989; Roux, 1999). In this
case, the electrical potential energy of the pore contents, ,
depends only on its net charge and dipole moment (Schu-
maker et al., 2000b). The energy is given by

qVI/2E , (23)
where q is the net charge of the pore contents,  is the
dipole moment calculated about z  0, and E  VI/L is the
constant electric field. For simplicity, we ignore the feath-
ering of the electric field that occurs at the channel en-
trances in the calculations of Jordan and Roux.
Proton transport through the channel corresponds to an
elementary charge moving through the thickness of the
membrane. This combination has units of dipole moment
and can be associated with dipole moment changes due to
diffusion across the proton and defect segments of the state
diagram in Fig. 1 D. The electrical potential energy drop
corresponding to ion or defect translocation through the
pore is proportional to the change in pore dipole moment.
This total change is 2A
H

15.86 e0Å in the case of proton
translocation and 2A
d

7.04 e0Å in the case of defect
translocation. Thus, 69% of the applied potential drop
drives proton translocation, and the rest drives the defects.
Framework model
Schumaker et al. (2000b) have constructed a framework
model for single-proton conduction through gramicidin
Single H Conduction through Gramicidin 2849
Biophysical Journal 79(6) 2840–2857

based on the mechanism of Fig. 1 C. The word “framework”
refers to the model’s design, which incorporates the results
of the molecular dynamics simulations. Mathematically, it
has the form of a pair of Nernst-Planck equations coupled
through their boundary conditions. Nernst-Planck equations
are the first integrals of Smoluchowski equations, which
model diffusion in the presence of systematic forces.
A state diagram for the model is presented in Fig. 1 D.
The top horizontal line segment represents the possible
proton states corresponding to the top row of Fig. 1 C. The
bottom horizontal line segment corresponds to the possible
defect states corresponding to the bottom row of Fig. 1 C.
Proton entrance and exit were not simulated by the molec-
ular dynamics. The single-proton model incorporates a very
simple representation of these processes.
Lumped state approximation of boundary regions
Define boundary regions on either side of the defect seg-
ment as regions receptive to protons entering from the
appropriate side. For example, protons may enter the chan-
nel from side I in Fig. 1 D if the defect dipole moment is in
boundary region I of Fig. 2 B. Mathematically, the boundary
regions are lumped as endpoints of the defect segment of the
state diagram. These endpoints are the boundary states bI
and bII of Fig. 1 D. Proton entrance occurs directly from the
boundary states. Schumaker et al. (2000b) show how
boundary conditions corresponding to the lumped state ap-
proximation are constructed.
Consider the state diagram of Fig. 1 D. Proton entrances
into the pore from side I occur with rate QI
b
ICI, where QI
b
is the probability that the boundary state bI is occupied and
ICI is the entrance rate of protons on side I. CI is the
concentration of the excess protons on side I, and I is a
second-order rate constant. A similar expression holds on
side II. The rates of transitions from the proton interval to
the boundary states are proportional to the proton probabil-
ity density at the endpoints of the proton interval. All
transitions satisfy detailed balance at equilibrium.
Through the boundary conditions on the Nernst-Planck
equations, it is possible to adjust the weight of the boundary
states bI and bII to be proportional to an integral of the
Boltzmann factor, exp [d(d)/kBT], over the correspond-
ing boundary regions of the d axis. Here, d is the defect
potential of mean force shown in Fig. 2 B. As a result, the
boundary states give the appropriate description of thermo-
dynamic equilibrium when the transmembrane electrical
potential difference is zero.
One limitation of the lumped state approximation of the
boundary regions is that the boundary states must act elec-
trically like points. Their electrical coordinates are adjust-
able parameters of the theory and are chosen to give the best
comparison with the exact mean first passage times over the
defect barrier. In Methods we compare estimates of the
mean first passage time for diffusion over the defect poten-
tial barrier as a function of applied potential. This analysis
shows that the lumped state approximation accurately mod-
els the dependence of mean first passage time on applied
potential.
Single-proton trajectories
The framework model is constructed as the limit of a
sequence of random walks (Schumaker et al., 2000b). The
sequence of random walks is indexed by an integer n,
denoting the number of states into which the proton and
defect segments are divided. Transition probabilities be-
tween the states are a function of n, and the framework
model is obtained in the limit n3 . Numerical simulations
suggest that transport described by the random walks accu-
rately simulates the limiting diffusion process for moderate
values of n. Here we use the random walks for n  57 to
illustrate ion trajectories underlying the framework model.
In particular, we illustrate the boundary behavior of their
trajectories.
In Fig. 6 A a short 50-ps segment of a proton trajectory is
shown. Proton escape is not limited by the shallow potential
energy well of the proton potential of mean force, but rather
by the boundary conditions. Fig. 6 B shows a trajectory on
the defect segment of the framework model state diagram
shown in Fig. 1 D. The extreme values of d correspond to
the boundary states at C
d . The trajectory makes rapid
excursions into the interior of the defect segment. Fig. 6 C
shows trajectories diffusing over the defect potential profile
on its entire interval of definition. The lumped state approx-
imation is not used. Qualitatively, the appearance of Fig. 6,
B and C, is somewhat similar in the region d  C
d . The
mean first passage times for transport over the defect barrier
are quantitatively very similar, as shown in Fig. 5 and
discussed in Methods.
Comparison with data of Eisenman et al.
Our main purpose in constructing the framework model of
single-proton conduction is to compare the results of the
molecular dynamics simulations (Pome`s and Roux, manu-
script in preparation) with experiment. Eisenman et al.
(1980) summarize several studies of ion conductance
through gramicidin, including an extensive investigation of
proton conductance for pH 1.7. Decker and Levitt (1988),
Akeson and Deamer (1991), Phillips et al. (1999), and
recently Cukierman (2000) have also made single-channel
measurements in this concentration range.
Eisenman et al. provide graphs of the logarithm of con-
ductance, log10 G50, and conductance ratio, G100/G0, as a
function of concentration C, where GV(C) is the proton
conductance of gramicidin A measured at symmetrical con-
centration C and applied potential V. We show these data in
Fig. 7, A and B. All of the data collected for pH  1.5 are
shown except for one outlier on each plot. To reduce clutter,
2850 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

several data points with pH  1.5 are omitted. The data
shown here fall on or close to smooth curves drawn through
the observations in the original paper.
Also shown in Fig. 7, A and B, are three fits of the
framework model to the data for pH  1.7. The data include
five conductances and eight conductance ratios. Each fit is
obtained by optimizing the parameters ta and  as described
in Methods; these parameters control the entrance and exit
rate of protons. The fits given by the solid and short-dashed
curves use the molecular dynamics estimates of the proton
and defect potentials of mean force. For the fit given by the
long-dashed curves, the defect potential of mean force
shown by Fig. 2 B was scaled by a constant factor to
increase the amplitude of the central barrier by 2kBT. All
three fits use the molecular dynamics estimates of the pro-
ton and defect diffusion coefficients, and a value of w 
2.4 Debye for the dipole moment of water, consistent with
the results of Duca and Jordan (1997).
The fits given by the solid and long-dashed curves use the
boundary regions shown in Fig. 2 B. We shall refer to these
as the wide boundary regions. The fit given by the short-
dashed curves uses boundary regions that are half as wide as
those shown in Fig. 2 B. We shall refer to these as narrow
boundary regions. Values of A
d , B
d , and C
d are given in the
legend of Fig. 2.
Conductance
Eisenman et al. (1980) measured single-channel proton con-
ductances in symmetrical solutions at 50 mV applied po-
tential. These data are shown as the filled circles in Fig. 7 A.
At the lowest concentrations, conductance is proportional to
concentration, corresponding to a slope of 1 on the log-log
plot. This is consistent with the idea that the proton entrance
is rate limiting, which is also strongly supported by Decker
and Levitt’s study (1988) of the effects of weak acids on
conductance at pH 2.75 and pH 3.75.
A striking feature of these data is the well-defined shoul-
der near pH 1.7. The existence of this shoulder has been
confirmed by Akeson and Deamer (1991), and very recently
a similar well-defined shoulder has been found in conduc-
tances of the RR dioxolane-linked gramicidin (Schumaker
et al., 2000a; Cukierman, 2000). A formula for the frame-
work model conductance is obtained by Schumaker et al.
(2000b); its form, for the case of symmetrical concentra-
tions and VI  0, is given by Eq. 21. Independent of the
details of the potentials of mean force or diffusion coeffi-
cients, one finds that there can be no shoulder such as that
displayed by the data in Fig. 7 A. One possible interpretation
of the experimental shoulder is that it signifies the onset of
significant two-proton conduction.
The quadratic dependence of the denominator of Eq. 21
on the symmetrical concentration C is unusual for a single-
ion conduction model. At low concentrations the conduc-
tance is similar to a Langmuir isotherm; it first increases
linearly with concentration and then begins to saturate. But
instead of asymptotically approaching a maximum value at
higher concentrations, the conductance given by Eq. 21
eventually decreases.
Physically, the decrease in conductance is due to compe-
tition between the time scales associated with defect diffu-
sion over the central barrier of Fig. 2 B and proton entrance.
Once a proton exits the channel in Fig. 1 D, a defect must
diffuse over the central barrier to complete the cycle around
the diagram associated with the net permeation of one
proton. If a proton first reenters from the same side as the
most recent exit, this cycling is frustrated. This is similar to
the clogging mechanism described by Schumaker and
MacKinnon (1990).
Each of the model curves shown in Fig. 7 A interpolates
the linearly increasing conductance at high pH very well,
but a difference arises in the region of the shoulder. The
model curve obtained by increasing the amplitude of the
defect segment central barrier is below the other two for
pH  2. This curve fits the data for pH  2 well, but if the
amplitude of the defect central barrier is increased much
further, it is no longer possible to interpolate the data near
the shoulder. This is shown by the sensitivity analysis given
below. A similar comparison between model curves and the
data arises if we use the defect potential of mean force
obtained from molecular dynamics but decrease the defect
diffusion coefficient. A good fit to the data is obtained when
the defect diffusion coefficient is a factor of 4 smaller than
that given by the molecular dynamics, but it becomes im-
possible to interpolate the shoulder if the diffusion coeffi-
cient is decreased much further.
The significance of the fits near the shoulder is not
completely clear because the single-proton model does not
account for the region above the shoulder. For example, if
the experimental shoulder is in fact due to the onset of
two-proton conductance, then this mechanism may make
some contribution below the shoulder as well, but this is not
taken into account here.
Conductance ratio
Eisenman et al. also measured the ratio of the conductance
at 100 mV to the conductance estimated near 0 mV. These
conductance ratios are shown as the filled circles in Fig. 7
B. A ratio greater than 1 corresponds to a superlinear I–V
curve, and a ratio less than 1 corresponds to a sublinear I–V
curve. These data have three very interesting features: a
limiting sublinear conductance ratio at high pH, an increas-
ing conductance ratio as pH decreases to 1, and then a
maximum ratio near pH 1, just above the shoulder of Fig. 7
A. Consistently, Cukierman has found a maximum in the
conductance ratio just above the shoulder in the RR con-
ductance data (Schumaker et al., 2000a).
The limiting ratio in the data of Eisenman et al. is
significantly less than 1, corresponding to sublinear I–V
Single H Conduction through Gramicidin 2851
Biophysical Journal 79(6) 2840–2857

curves. This suggests that the rate of proton entrance is not
strongly dependent on the applied potential VI. As a sim-
plifying assumption, our framework model assumes that
entrance is independent of VI (Schumaker et al., 2000b).
Indeed, the model I–V curves are also sublinear at high pH.
There is, however, a difference in limiting conductance
ratios between the solid and long-dashed curves, which
were constructed using wide boundary regions, and the
short-dashed curve, which was constructed using narrow
boundary regions. The limiting ratio is larger for the narrow
boundary regions. Similarly, the limiting conductance ratio
depends on the assumed value w of the dipole moment of
a water molecule. The limiting conductance ratio increases
as w increases. This interplay between the width of the
boundary regions and the value of w is also revealed by the
sensitivity analysis, which we discuss next. Each of the
model fits interpolates the data well in the regime of in-
creasing conductance ratios.
Sensitivity analysis
In this section we test the sensitivity of the optimized fits of
the single-proton model to the data of Eisenman et al. as we
perturb the height of the defect barrier and the value of the
defect diffusion coefficient obtained from molecular dy-
namics, and as we also perturb the assumed value of the
water dipole moment, w. This analysis is motivated in two
ways. First, the data we are fitting are limited: just five
conductance data points and eight curve shape data points
for pH  1.7. To what extent does optimization of ta and 
give the model a general capability to interpolate these data
points? Can we arbitrarily vary the potential profile and
diffusion coefficients and obtain equally good fits? Second,
the molecular dynamics simulations did not include a rep-
resentation of the membrane or the aqueous solution on
either side. Therefore, the potentials of mean force we use
neglect contributions due to interactions with these compo-
nents of the channel environment. If these interactions
change the height of the defect potential barrier, shown in
Fig. 2 B, this may conceivably be reflected in an improved
optimized fit with the data when the height is perturbed.
To test the sensitivity of the optimized fit to the height of
the defect barrier, we ran a series of analyses, using modi-
fied defect potentials of mean force. The error (the sum of
the squared differences between the conductance data points
and the model plus twice the sum of the squared differences
between the curve shape data points and the model) was
computed as a function of the change in the height of the
defect potential barrier. The defect potential was rescaled to
change the amplitude of the central barrier by increments of
2kBT. This was done for both the wide boundary regions
shown in Fig. 2 B and for narrow boundary regions de-
scribed in the legend of Fig. 2. Results are shown in Fig. 8
A. The height of the defect barrier can be greatly decreased,
and equally good fits are obtained. However, if this height
FIGURE 8 Sensitivity analysis. The error is a measure of the deviation
between the model and the Eisenman conductance data for pH  1.7 and
between the model and the Eisenman curve shape data for pH  1.7. It is
scaled so that the fit using the molecular dynamics potential of mean force
and diffusion coefficients with wide boundary regions and w  2.4 Debye
has error  1. Diamonds connected by solid lines show the error calculated
using wide boundary regions, and stars connected by dotted lines show the
error calculated using narrow boundary regions. (A) Error as a function of
change in barrier height. The amplitude of the molecular dynamics poten-
tial of mean force was rescaled to several different barrier heights. The
abscissa is the difference between the rescaled barrier height and the
original molecular dynamics barrier height; 0 corresponds to the molecular
dynamics potential. (B) Error as a function of defect diffusion coefficient.
The diffusion coefficient obtained from molecular dynamics is multiplied
by a prefactor p. The abscissa is log2 p; 0 corresponds to the molecular
dynamics diffusion coefficient. (C) Error as a function of the assumed
dipole moment of water.
2852 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

is increased by much more than 2kBT, good fits to the data
are no longer possible.
To test the sensitivity of the optimized fit to the values of
the defect diffusion coefficient, we also ran a series of
analyses using modified diffusion coefficients. These may
be expressed in the form pd, where d is the result from
molecular dynamics and p is a prefactor. Fig. 8 B shows the
same measure of error plotted as a function of log2 p.
Results for both wide and narrow boundary regions are
shown. The defect diffusion coefficient can be greatly in-
creased beyond the value obtained from molecular dynam-
ics, and equally good fits are obtained. However, if this
diffusion coefficient is decreased by much more than a
factor of 4, the fits rapidly deteriorate.
The parameter ta is inversely proportional to the proton
exit rate. When either the defect diffusion coefficient is
increased or the barrier height is decreased from the molec-
ular dynamics values, the optimized values of ta vary little.
These results suggest that, for these series, transport across
the defect potential of mean force barrier is not rate limiting.
However, when either the defect diffusion coefficient is
decreased or the barrier height is increased, the optimized
values of ta rapidly decrease, corresponding to a rapidly
increasing optimal exit rate. Thus, for these series, the rates
of proton exit and transport across the defect barrier are both
important.
The proton potential of mean force, H, obtained from
the molecular dynamics simulations, is a shallow well. The
proton diffusion coefficient, DH, obtained from the Laplace
transform analysis of the simulated velocity autocorrelation
functions, is very high. Consequently, proton translocation
through the channel is not rate limiting. Our results are
insensitive to moderate changes in the proton potential
profile and diffusion coefficient.
Fig. 8 C shows how the optimized fit varies with the
assumed value of the dipole moment of water. Fits to the
Eisenman et al. (1980) data set using the narrow boundary
regions are best near w  2.4 Debye, whereas fits using
the wide boundary regions are best for w
3 Debye. This
result is somewhat surprising. The wide boundary regions,
which contain the potential minima at d 6.5 e0Å, seem
more reasonable than the narrow boundary regions, which
are intervals of w beyond the minima. However, Duca and
Jordan only found values of w as high as 3 Debye for
waters that are nearest neighbors of a cation in the channel.
A close analysis of the fits shown in Fig. 7 B suggests that
the largest difference between the fits lies in the limiting
conductance ratios at low concentrations. The conductance
ratios for the wide boundary regions lie below most of the
data points. This suggests a speculative explanation for the
surprising result. Interfacial polarization would be expected
to increase the conductance of the experimental system at
low concentrations and high applied potentials (Everitt and
Haydon, 1968; Walz et al., 1969), an effect that is not taken
into account in our model. The biggest effect on the Eisen-
man data set would be to increase values of G100 at low
concentrations, increasing the experimental conductance ra-
tio in that regime.
Comparison with other data
Decker and Levitt (1988) have studied proton conductance
through gramicidin, publishing I–V curves at pH 3.75 and
2.75. Fig. 7, C and D, compares their data to the fits of the
three single-proton models considered in Fig. 7, A and B.
The good agreement between the data and the model curves
suggests that the results of Decker and Levitt are consistent
with those of Eisenman et al. (1980). The data do show
somewhat larger currents than the model, particularly at
large VI. These may well be due to small differences in
experimental conditions. But the difference is also consis-
tent with increased conductance in the experimental system
due to interfacial polarization, which leads to I–V curves
that are asymptotic to lines of positive slope at high poten-
tials. This effect is very pronounced for Cs conductance
through gramicidin at 10 mM (Andersen, 1983).
Decker and Levitt (1988) also calculate a capture radius,
rc, for protons based on their data. The capture radius is
defined by the formula for the flux due to diffusion toward
an absorbing hemisphere of radius rc. Physically this flux
would correspond to a saturating current through a channel
at low bulk concentration (where the current is entrance-
limited) and at high applied potential (where the boundary
condition at the high potential side of the membrane could
be reasonably approximated as perfectly absorbing). Based
on their current at 200 mV applied potential and at pH 3.75
(Fig. 7 C), they estimated rc  0.87 Å.
We have estimated rc by using the single-proton model
formula for the saturating current in the limit of high applied
potentials (Schumaker et al., 2000b). For the molecular
dynamics parameter set and wide boundary regions (Fig. 7,
solid curves) we obtain rc  0.90 Å. When the amplitude of
the central barrier is increased by 2kBT (long-dashed
curves), we also obtain rc  0.90 Å. For narrow boundary
regions (short-dashed curves) we obtain rc  2.29 Å. The
higher values of rc for narrow boundary regions reflect the
fact that, to fit the data, the single-proton model entrance
rate must be increased to compensate as the boundary
regions are made more narrow. This leads to a higher
saturating current at very high applied potentials.
Akeson and Deamer (1991) have also reported single-
channel measurements of proton conduction through gram-
icidin. Their measurements of proton currents at 59 mV as
a function of pH are in qualitative agreement with the
findings of Eisenman et al., including a shoulder in conduc-
tance between pH 1 and pH 2. However, they also publish
an I–V curve at pH 2, which shows complete saturation of
current as a function of voltage above VI  125 mV. The
shapes of our corresponding I–V curves are qualitatively
similar to those shown in Fig. 7, C and D; they do not have
Single H Conduction through Gramicidin 2853
Biophysical Journal 79(6) 2840–2857

the sharply saturating character of the data reported by
Akeson and Deamer.
SUMMARY
Single-proton model
We summarize the construction of the single-proton con-
duction model. Pome`s and Roux (1997) performed two sets
of molecular dynamics simulations of waters in the grami-
cidin pore. In one set the pore includes an excess proton,
and we refer to the resulting proton potential of mean force.
This is shown in Fig. 2 A as a function of the orientation
moment of the pore contents. It has the form of a shallow
well. The molecular dynamics also estimates a diffusion
coefficient for proton transport in terms of the dipole mo-
ment reaction coordinate. It is 40 times the diffusion
coefficient of protons in water, or about twice the diffusion
coefficient of protons in ice.
In the second set of simulations, there is no excess proton.
It is energetically favorable for the pore waters to be aligned
with similar components of the dipole moment parallel to
the pore axis. The alignment of the pore waters reverses
when a packing defect propagates through the water chain.
The corresponding defect potential of mean force, as a
function of pore dipole moment, is shown in Fig. 2 B. It has
a central potential barrier, 6kBT high, surrounded by min-
ima. An effective diffusion coefficient for defects is also
estimated by the molecular dynamics. Its magnitude is con-
sistent with the rate of molecular tumbling in bulk water.
The molecular dynamics simulations described above did
not model proton entrance into the pore. To use these results
in a theory of proton permeation through gramicidin, it is
necessary to add a mechanism for proton entrance. Schu-
maker et al. (2000b) have developed a framework model of
proton conduction through gramicidin. This probabilistic
model is designed to incorporate the results from molecular
dynamics, and includes a simple description of ion entrance
and exit.
The physical assumption made by the framework model
is that proton entrance is possible only when the dipole
moment of the pore water column favors proton entry. This
idea is illustrated by the cartoons of Fig. 1 C. A state
diagram for the mechanism is shown in Fig. 1 D. The top
segment of the diagram represents the set of possible proton
occupied states. The bottom segment includes a continuum
of states representing defect transport over the central bar-
rier in Fig. 2 B and endpoints at both ends. The endpoints
represent defect occupation of boundary regions near the
potential minima.
The endpoints are constructed so that, in a state of sym-
metrical equilibrium, their probability is equal to the prob-
ability of defect occupation of the boundary regions (Schu-
maker et al., 2000b). We call this hybrid of discrete and
continuous representations the lumped state approximation.
The endpoints, which are the lumped states, must respond to
an applied electrical field like points. As long as the bound-
ary regions are not too wide, their effective electrical coor-
dinate can be adjusted so that the mean first passage times
to cross the central barrier using the lumped states agree
closely with estimates that do not involve this approxima-
tion. This is demonstrated by the calculations used to con-
struct Fig. 5.
The lumped state approximation makes it possible to
solve the framework model analytically (Schumaker et al.,
2000b). Analytical expressions for the current and conduc-
tance are found. These include dependence on the proton
and defect potentials of mean force as arbitrary functions.
Our model for single-proton conduction through gramicidin
is the framework model combined with the results from
molecular dynamics. It has three free parameters, control-
ling the rate of ion entrance, ion exit, and the width of the
boundary regions.
Comparison with experiment
Eisenman et al. (1980) summarized an extensive study of
proton conductance in gramicidin. The data in Fig. 7, A and
B, are taken from similar figures in their paper. The con-
ductance data in Fig. 7 A have a well-defined shoulder for
1  pH  2, with an inflection point near pH 1.7. This
shoulder is impossible to fit using the single-proton model,
as can be seen from the analytical formula for conductance
(Eq. 21). The existence of the shoulder has been confirmed
in gramicidin A by Akeson and Deamer, and a similar
feature has been seen recently in the RR dioxolane dimer
(Cukierman, 2000). Furthermore, a maximum in the RR
curve shape ratios near [H]  100 mM has been observed
that is qualitatively similar to that seen in the data of
Eisenman et al., shown in Fig. 7 B.
Our interpretation is that the inflection point signifies the
onset of two-proton conductance, in agreement with Eisen-
man et al. We therefore fit the single-proton conduction
model to the data in Fig. 7, A and B, only for pH  1.7. To
make the fits, we consider two values for the width of the
boundary regions. The wide boundary regions are shown in
Fig. 2 B, while the narrow boundary regions have only half
the width shown. For a fixed choice of boundary regions,
the two parameters controlling the proton entrance and exit
rate are optimized. Good fits to the data are achieved, as
shown in Fig. 7, A and B. We have studied their significance
by perturbing the defect potential of mean force and diffu-
sion coefficient. Fig. 8 A was constructed by repeating the
optimization procedure using rescaled defect potentials with
a range of barrier heights. The barrier height obtained from
molecular dynamics is within a range of values that give
nearly optimal fits. A similar result was found by varying
the value of the defect diffusion coefficient, as shown by
Fig. 8 B.
2854 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

An important limitation of the analysis given here is that
only single-proton transport through gramicidin is de-
scribed, and we have no quantitative comparison with the
data for pH  1.7. However, we are currently extending the
present model by adding rate-theory-like terms describing
transitions into and out of a doubly occupied state. A good
fit to the Eisenman et al. conductance data on both sides of
the shoulder shown in Fig. 7 A can be achieved (Schumaker
et al., 2000a). DeCoursey and Cherny (1999) have argued
that multiple proton occupancy of the gramicidin pore is
unlikely because it would require an energetically expensive
defect in the hydrogen-bonded chain between the two cen-
ters of excess charge. Indeed, the shoulder in Fig. 7 A,
between two regimes in which conductance is proportional
to [H], may be interpreted as corresponding to an in-
creased energy barrier for ion entrance into the channel.
Molecular dynamics studies would provide useful insights
into the energetics of double proton occupancy and the
relationship between the second proton entrance/exit pro-
cess and the configuration of the hydrogen bonded chain.
Phillips et al. (1999) have recently studied proton con-
ductance in side-chain analogs of gramicidin A channels.
Fluorination of the channel Trp11, Trp13, or Trp15 side
chains is found to inhibit proton transport, and replacement
of one or more Trps with Phe enhances transport. The
observations with the fluorinated analogs were made at
three proton concentrations in the range 0.018 N [H]
0.18 N, and those with the Phe analogs were made at
[H]  0.1 N. They can be interpreted consistently as
suggesting that increasing the electrical potential of the
channel interior with respect to the membrane surface has
the effect of increasing conductance. Phillips et al. explain
these results by proposing that the rate-limiting step may be
water reorientation mediated by an exit-initiated defect car-
rying a partial negative charge across the pore interior. This
seems natural in view of the fact that Pome`s and Roux
(1997, and manuscript in preparation) found a substantial
energy barrier to water reorientation (shown in Fig. 2 B).
Increasing the electrical potential in the pore interior would
lower this barrier for a negatively charged defect. Phillips et
al. further argue that the dipole-dipole interaction between
tryptophan and the pore waters stabilizes the orientation of
waters at the side of the channel from which a proton has
just exited. Then the Trp3 Phe substitution should increase
conductance if the creation of defects in the hydrogen-
bonded chain of waters at the exit side were a rate-limiting
step.
However, the simulations of Pome`s and Roux (1996;
manuscript in preparation) suggest that waters at either end
of the pore are somewhat disordered. This disorder is asso-
ciated with the minima at either end of the defect potential
of mean force in Fig. 2 B. These results suggest that defect
creation at the ends of the water chain is not rate limiting.
Furthermore, there may be a difficulty with reorienting the
channel waters by beginning at the end of the pore from
which the proton has most recently exited. Between the two
subcolumns of aligned water molecules there would be
either a Bjerrum L defect with two water oxygens opposed,
or a similar configuration with a water interposed between
the subcolumns and making a hydrogen bond with each.
Apparently, one or two of the waters in these conformations
participate in only two hydrogen bonds (novel modes of
hydrogen bonding may conceivably alter this conclusion;
see Vargas et al., 2000). In contrast, if water reorientation
begins from the other side of the pore, the defect can
propagate through a series of conformations in which every
water makes at least three hydrogen bonds (Pome`s, 1999;
Pome`s and Roux, manuscript in preparation). Thus there
may be an energy penalty associated with the water reori-
entation mechanism envisioned by Phillips et al.
The present model provides an alternative interpretation
of the experimental findings of Phillips et al. In the region
of the shoulder of the Eisenman conduction data, our results
suggest that proton exit is a rate-limiting step. Specifically,
the sensitivity analysis shows that proton exit and water
reorientation rates are comparable in the shoulder region if
water reorientation is somewhat slower than that calculated
by the molecular dynamics (fits become impossible if water
reorientation is much slower). The proton exit rate is lim-
iting by itself in the shoulder region if water reorientation is
much faster than that calculated by the molecular dynamics.
Therefore, an increased electrical potential in the channel
interior might plausibly speed proton exit and thus increase
conductance, which is consistent with the observations. It
would be very interesting to study the conductance of the
side-chain analogs for a wider range of concentrations.
APPENDIX A: ESTIMATE OF DIPOLE DIFFUSION
COEFFICIENT FROM ROTATIONAL TUMBLING
A very simple estimate of the rotational diffusion coefficient of water may
be developed from considerations independent of the molecular dynamics
calculation. Consider Brownian motion of an overdamped degree of free-
dom  with friction coefficient  in a harmonic potential of spring constant
K. By equipartition, the mean value of the coordinate is given by 2 
kBT/K. The time correlation function is
t
0

2exp t/
, (24)
where
 /K is the time constant. Its time derivative is
d/dtt
0
t0
 kBT/
i , (25)
where the Einstein relationship between friction and diffusion coefficient

i is used. The pore of gramicidin contains 10 water molecules. The total
dipole moment is



m1
10
m , (26)
and the time correlation function is
t
0



m1
10

n1
10
mt
n0
 . (27)
Single H Conduction through Gramicidin 2855
Biophysical Journal 79(6) 2840–2857

For the purpose of making a simple estimate, assume complete indepen-
dence of the water dipoles,
mt
n0

mnmt
n0
 , (28)
where mn  1 if m  n and mn  0 otherwise. Insert this expression into
Eq. 27 to obtain
t
0

102et/
. (29)
Using Eq. 25, we then get

i

102/
. (30)
The dipole moment of a water molecule in a gramicidin-like pore is
estimated to be 0.50 e0Å (Duca and Jordan, 1997), and the rotational
tumbling relaxation time in bulk solution is 7 ps at 20°C (Eisenberg and
Kauzmann, 1969), giving the estimate

i
 0.36e0 Å
2ps
1 . (31)
APPENDIX B: MFPT TO CROSS DEFECT
BARRIER
In this appendix we obtain two expressions for the mean time for the defect
to cross the central barrier shown in Fig. 2 B. Neither of these expressions
assumes the lumped state approximation used by the single-proton con-
duction model.
Exact mean first passage time
Calculation of mean first passage time for a diffusion process (represented
by a Smoluchowski or Nernst-Planck equation) to escape a region with
given boundary conditions is a standard topic in the theory of stochastic
processes (Karlin and Taylor, 1981). A derivation based on the theory of
random walks is given in the appendix of Schumaker and Kentler (1998).
A related textbook discussion is given by Zauderer (1989).
Consider diffusion over the total defect potential W()  d() 

d(), where d is the intrinsic component of the potential, pictured in
Fig. 2 B, and d()  dE is the component due to the applied
transmembrane potential. Let t() be the mean time for a diffuser begin-
ning at reaction coordinate   dd to first reach the potential minimum
at min
d
 6.5 e0Å, given a reflecting boundary at A
d
 8.1 e0Å. The
mean first passage time satisfies
D


2t

2 
W
kBT
t



 1 , (32)
where D is the associated diffusion coefficient and the boundary conditions
are t (A
d )  0 and t(min
d )  0, where min
d
 
d
min
d . The solution may
be found by the method of variation of parameters. It is
Dt


A
d
min
d
smin
d
 s
eW
/kBTd


A
d

s
 s
eW
/kBTd, (33)
with s given by
s


A
d

eW
/kBTd . (34)
Then t(min
d ) is the mean first passage time to first reach the potential
minimum at d  min
d , given that the trajectory starts at the potential
minimum at d  min
d and that there is a reflecting boundary at d 
A
d . This is the exact mean first passage time referred to in the Methods.
Kramers’ approximation
An approximate, but simpler, expression for the rate of barrier crossing is
given by Kramers’ escape theory (e.g., Risken, 1989). This theory assumes
that the diffuser attains a quasiequilibrium state in a deep potential well
before making a transit over a high barrier. Under these conditions, transits
will be nearly exponentially distributed in time. The transition rate from
Kramers’ theory is
k
D/


3
4
eW
/kBTd

1
2
eW
/kBTd


, (35)
and the corresponding mean first passage time is 1/k. In this expression, the
pair 1 and 2 span the potential minimum at 3, and the pair 3 and 4
span the potential maximum. The process is absorbed at 4. To construct
Fig. 5, these integrals are evaluated numerically. When the top of the
barrier and the bottom of the well can be approximated by quadratic
extrema, the classical escape rate theory expression is obtained (e.g.,
Cooper et al., 1988; Risken, 1989).
MFS enjoyed conversations and communications with Mark Akeson, See-
Wing Chiu, Sam Cukierman, Tom DeCoursey, and Eric Jakobsson, as well
as the hospitality of the Departments of Chemistry and Physics at the
University of Montreal, where this work began. He gratefully acknowl-
edges the help of both Marc Souille with the computer system at the
University of Montreal and of Wonpil Im with CHARMM calculations.
MFS also enjoyed the hospitality extended by Brigham Young University,
and particularly thanks Dr. David Busath for extensive and enjoyable
conversations. He finally thanks Ann Schumaker for carefully editing the
manuscript.
MFS is supported by grant MCB 9630475 from the National Science
Foundation. RP is supported in part by the U.S. Department of Energy
through the Los Alamos National Laboratory Directed Research and De-
velopment Grant for Bioremediation. BR is supported by the Medical
Research Council of Canada.
REFERENCES
Akeson, M., and D. W. Deamer. 1991. Proton conduction by the gramici-
din water wire. Biophys. J. 60:101–109.
Andersen, O. S. 1983. Ion movement through gramicidin A channels.
Interfacial polarization effects on single-channel current measurements.
Biophys. J. 41:135–146.
Arseniev, A. S., I. L. Barsukov, V. F. Bystrov, A. L. Lomize, and Yu. A.
Ovchinnikov. 1985. 1H-NMR study of gramicidin-A transmembrane ion
channel: head-to-head right handed single stranded helices. FEBS Lett.
186:168–174.
Berne, B. J., M. Borkovec, and J. E. Sraub. 1988. Classical and modern
methods in reaction rate theory. J. Phys. Chem. 92:3711–3725.
Berne, B. J., and R. Pecora. 1976. Dynamic Light Scattering. Wiley, New
York.
2856 Schumaker et al.
Biophysical Journal 79(6) 2840–2857

Boyer, P. 1997. The ATP synthase—a splendid molecular machine. Annu.
Rev. Biochem. 66:717–749.
Busath, D. D., C. D. Thulin, R. W. Hendershot, L. R. Phillips, P. Maughan,
C. D. Cole, N. C. Bingham, S. Morrison, L. C. Baird, R. J. Hendershot,
M. Cotten, and T. A. Cross. 1998. Noncontact dipole effects on channel
permeation. I. Experiments with (5F-indole)Trp13 gramicidin A chan-
nels. Biophys. J. 75:2830–2844.
Chang, G., R. H. Spencer, A. T. Lee, M. T. Barclay, and D. C. Rees. 1998.
Structure of the MscL homolog from Mycobacterium tuberculosis: a
gated mechanosensitive ion channel. Science. 282:2220–2230.
Cooper, K. E., P. Y. Gates, and R. S. Eisenberg. 1988. Diffusion theory and
discrete rate constants in ion permeation. J. Membr. Biol. 106:95–105.
Cooper, K., E. Jakobsson, and P. Wolynes. 1985. The theory of ion
transport through membrane channels. Prog. Biophys. Mol. Biol. 46:
51–96.
Cowan, S. W., T. Schirmer, G. Rummel, M. Steiert, R. Ghosh, R. A.
Pauptit, J. N. Jansonius, and J. P. Rosenbusch. 1992. Crystal structures
explain functional properties of two E. coli porins. Nature. 358:727–733.
Crouzy, S., T. B. Woolf, and B. Roux. 1994. A molecular dynamics study
of gating in dioxolane linked gramicidin-A channels. Biophys. J. 67:
1370–1386.
Cukierman, S. 2000. Proton mobilities in water and in different stereoiso-
mers of covalently linked gramicidin-A channels. Biophys. J. 78:
1825–1834.
Decker, E. R., and D. G. Levitt. 1988. Use of weak acids to determine the
bulk diffusion limitation of H ion conductance through the gramicidin
channel. Biophys. J. 53:25–32.
DeCoursey, T. E., and V. V. Cherny. 1999. An electrophysiological com-
parison of voltage-gated proton channels, other ion channels, and other
proton channels. Isr. J. Chem. 39:409–418.
Doyle, D. A., J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L.
Cohen, B. T. Chait, and R. MacKinnon. 1998. The structure of the
potassium channel: molecular basis of K conduction and selectivity.
Science. 280:68–76.
Duca, K. A., and P. C. Jordan. 1997. Ion-water and water-water interac-
tions in a gramicidinlike channel: effects due to group polarizability and
backbone flexibility. Biophys. Chem. 65:123–141.
Eisenberg, D., and W. Kauzmann. 1969. The Structure and Properties of
Water. Oxford University Press, New York.
Eisenman, G., B. Enos, J. Ha¨gglund, and J. Sandbloom. 1980. Gramicidin
as an example of a single-filing ion channel. Ann. N.Y. Acad. Sci.
339:8–20.
Elston, T. C., and G. Oster. 1997. Protein turbines. I. The bacterial flagellar
motor. Biophys. J. 73:703–721.
Everitt, C. T., and P. A. Haydon. 1968. Electrical capacitance of a lipid
membrane separating two aqueous layers. J. Theor. Biol. 18:371–379.
Heinemann, S., and F. J. Sigworth. 1989. Estimation of Na dwell time in
the gramicidin A channel. Na ions as blockers as H currents. Biochim.
Biophys. Acta. 987:8–14.
Hille, B. 1992. Ionic Channels of Excitable Membranes, 2nd Ed. Sinauer
Associates, Sunderland, MA.
Jackson, J. D. 1975. Classical Electrodynamics. Wiley, New York.
Jordan, P. C. 1982. Electrostatic modeling of ion pores. Biophys. J.
39:157–164.
Jordan, P. C., R. J. Bacquet, J. A. McCammon, and P. Tran. 1989. How
electrolyte shielding influences the electrical potential in transmembrane
ion channels. Biophys. J. 55:1041–1052.
Karlin, S., and H. M. Taylor. 1981. A Second Course in Stochastic
Processes. Academic Press, New York.
Ketchem, R. R., W. Hu, and T. A. Cross. 1993. High-resolution confor-
mation of gramicidin-a in a lipid bilayer by solid state NMR. Science.
261:1457–1460.
Phillips, L. R., C. D. Cole, R. J. Hendershot, M. Cotten, T. A. Cross, and
D. D. Busath. 1999. Non-contact dipole effects on channel permeation.
III. Anomalous proton conductance effects in gramicidin. Biophys. J.
77:2492–2501.
Pome`s, R. 1999. Theoretical studies of the Grotthus mechanism in biolog-
ical proton wires. Israel J. Chem. 39:387–395.
Pome`s, R., and B. Roux. 1996. Structure and dynamics of a proton wire: a
theoretical study of H translocation along the single file water chain in
the gramicidin A channel. Biophys. J. 71:19–39.
Pome`s, R., and B. Roux. 1997. Free energy profiles governing H con-
duction in proton wires. Biophys. J. 72:A246.
Risken, H. 1989. The Fokker-Planck Equation. Springer-Verlag, New
York.
Roux, B. 1999. Statistical mechanical equilibrium theory of selective ion
channels. Biophys. J. 77:139–153.
Roux, B., and M. Karplus. 1993. Ion transport in the gramicidin channel:
free energy of the solvated right-hand dimer in a model membrane.
J. Am. Chem. Soc. 115:3250–3262.
Schumaker, M. F., and C. Kentler. 1998. Far-field analysis of coupled bulk
and boundary layer diffusion. Biophys. J. 74:2235–2248.
Schumaker, M. F., and R. MacKinnon. 1990. A simple model for multi-ion
permeation. Biophys. J. 58:975–984.
Schumaker, M. F., R. Pome`s, B. Roux, and S. Cukierman. 2000a. New
experimental results and an extended model of proton conduction
through gramicidin. Biophys. J. 78:348A.
Schumaker, M. F., R. Pome`s, and B. Roux. 2000b. A framework model for
single proton conduction through gramicidin. Biophys. J. In press.
Stillinger, F. H., and C. W. David. 1978. Polarization model for water and
its ion dissociation products. J. Phys. Chem. 69:1473–1484.
Straub, J., and B. J. Berne 1988. Molecular dynamics study of an isomer-
izing diatomic in a Lennard-Jones fluid. J. Phys. Chem. 89:4833–4847.
Vargas, R., J. Garza, D. Dixon, and P. H. Hay. 2000. How strong is the
C  H. . .0 | C hydrogen bond? J. Am. Chem. Soc. 122:4750–4755.
Walz, D., E. Bamberg, and P. La¨uger. 1969. Nonlinear electrical effects in
lipid bilayer membranes. I. Ion injection. Biophys. J. 9:1150–1159.
Weber, T. A., and F. H. Stillinger. 1982. Reactive collisions of H3SO4 and
OH studied with the polarization model. J. Chem. Phys. 86:
1314–1318.
Woolf, T. B., and B. Roux. 1994. The conformational flexibility of o-
phosphorylcholine and o-phosphiorylethanolamine: a molecular dynam-
ics study of solvation effects. J. Am. Chem. Soc. 116:5916–5926.
Zauderer, E. 1989. Partial Differential Equations of Applied Mathematics.
Wiley, New York.
Single H Conduction through Gramicidin 2857
Biophysical Journal 79(6) 2840–2857