Framework Model For 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, and ‡Groupe de Recherche en
Transport Membranaire, De´partements de Physique et de Chimie, Universite´ de Montre´al, Que´bec H3C 3J7, Canada
ABSTRACT This paper describes a framework model for proton conduction through gramicidin; a model designed to
incorporate information from molecular dynamics and use this to predict conductance properties. The state diagram
describes both motion of an excess proton within the pore as well as the reorientation of waters within the pore in the absence
of an excess proton. The model is constructed as the diffusion limit of a random walk, allowing control over the boundary
behavior of trajectories. Simple assumptions about the boundary behavior are made, which allow an analytical solution for the
proton current and conductance. This is compared with corresponding expressions from statistical mechanics. The random
walk construction allows diffusing trajectories underlying the model to be simulated in a simple way. Details of the numerical
algorithm are described.
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 numeral R
denotes either side I or side II of the channel.
Latin Symbols
a, aˆ Weights of boundary regions. Eqs.
38, 43, and 65.
bR Boundary state R. See Fig. 2 B.
CR, cR Bulk concentration on side R. See
Eqs. 32 and 57.
C0 A concentration introduced by Eq.
38; see also Eq. 72.
C• The unit concentration, e.g. 1M; see
Eq. 44.

s Diffusion coefficient of dipole
reaction coordinate. See above Eq.
18.
di Defect state i; see Fig. 2 C.
E Electric field in pore interior. See
above Eq. 3.
e0 Elementary electrical charge.
fH Proton electrical distance. See Fig.
2 A and Eq. 12.
f A
d , f B
d , f C
d
 f d Defect electrical distances. See Fig.
2 A and Eqs. 13–15.
gs Integral defined by Eq. 120.
Hi Proton state i; see Fig. 2 C.
hs Integral defined by Eq. 114.
is Integral defined by Eq. 121.
Js Flux of species s; see Eq. 31.
I Current through channel; see Eq. 83.
kB Boltzmann’s constant.
Ks Integration constants for Ps; see Eq.
34.
L Spatial length of channel. See the
first paragraph in Construction of the
Model.

s Length of species s reaction
coordinate interval. See above Eq.
18.
n Number of random walk gridpoints
for each species s. See Fig. 2 C.
Ps, ps Probability density at s. See Eq. 27
and above Eq. 92.
Qi
s Probability that state si is occupied;
see above Eq. 21.
QR
b Probability that bR is occupied; see
above Eq. 32.
QFW
H Framework model probability that a
proton occupies the channel; Eq. 68.
QSM
H Statistical mechanical probability of
proton occupation; see Eq. 67.
T Absolute temperature.
ta, ts Access time or characteristic time for
species s. See Eq. 55 and above.
VI Applied potential on side I; see
above Eq. 3.
Ws, ws Total energies of species s. See Eqs.
17 and above 92.
z Spatial coordinate coaxial with the
pore. See the first paragraph in
Construction of the Model.
Received for publication 14 February 2000 and in final form 2 October
2000.
Address reprint requests to Mark F. Schumaker, Washington State Uni-
versity, Department of Pure and Applied Mathematics, Pullman, Washing-
ton 99164. Tel.: 509-335-7273; Fax: 509-335-1188; E-mail: schumaker@
wsu.edu.
Re´gis Pome`s’ present address is Structural Biology and Biochemistry,
Hospital For Sick Children, Toronto, Ontario M5G 1X8, Canada.
Benoıˆt Roux’s present address is Department of Biochemistry and Struc-
tural Biology, Weill Medical College of Cornell University, 1300 York
Ave., New York, New York 10021.
© 2001 by the Biophysical Society
0006-3495/01/01/12/19 $2.00
12 Biophysical Journal Volume 80 January 2001 12–30

Greek Symbols
RCR Entrance transition probability on side R;
see above Eq. 32.
 (kBT)
1.
R Exit transition probability on side R; see
above Eq. 32.
i
s Forward transition probability from state i.
See Fig. 2 C and Eq. 18.
t Random walk time step; see Eq. 18.
 Factor in t independent of n, Eq. 20.
i
s Backward transition probability from state
i. See Fig. 2 C and Eq. 19.
 Defined by Eq. 44.
R Transition probability from bR to interior of
defect interval; see above Eq. 32.

s, i
s Dipole moment reaction coordinate for
species s. See Fig. 1 and below Eq. 17.
A
s Maximum extent of coordinate interval.
See Fig. 1.
B
d Effective electrical coordinates of boundary
regions. See Fig. 2 A and above Eq. 16.
C
d Maximum extent of interior of defect
interval. See above Eq. 1.

R Transition probability from interior of
defect interval to bR; see above Eq. 32.

s, i
s Dimensionless reaction coordinate for
species s. See Eqs. 90 and 91.
B
d See below Eq. 44.
C
d Eq. 2.
A
H Eq. 1.

s, s Potential of mean force for species s. See
Fig. 1 and above Eqs. 1 and 92.

s, s Relative potentials of mean force. See Eqs.
1 and above Eq. 92.

s,
s Applied electrostatic potential energy. See
Fig. 1 and Eq. 3 and above Eq. 92.
I,
I Energy of an elementary charge in potential
VI. See Eq. 37 and above Eq. 92.
INTRODUCTION
This paper describes the construction of a framework model
of single proton conduction through the ion channel gram-
icidin. A framework model is a kinetic model, designed to
incorporate potentials of mean force and diffusion coeffi-
cients computed by molecular dynamics simulations on a
very short time scale, and then use this information to
calculate conductances and associated observable quantities
measured on a much longer time scale. The reaction coor-
dinates of molecular dynamics simulations parameterize a
simplified configuration space for the system being mod-
eled. They explicitly parameterize the degrees of freedom
thought to be most important for describing the system. But
they average over fast variables; for example, those that
describe intermolecular vibrations. A framework model that
incorporates information from molecular dynamics is in the
same sense a simplified model of configuration space.
Framework models are somewhat similar to rate theory
models (reviewed by Hille, 1992), which can be regarded as
zero dimensional approximations of configuration space.
Rate theory models consist of states (e.g., the empty chan-
nel, an ion occupying one binding site, an ion occupying a
second binding site, etc.) and transitions between them. As
models of configuration space, they naturally incorporate
restrictions on internal degrees of freedom due to the nature
of condensed phase motion on molecular scales. For exam-
ple, one may easily describe channels whose occupancy is
limited to a single ion (e.g., La¨uger, 1973), or multiply
occupied channels (e.g., Hille and Schwarz, 1978). In this
respect, rate theory models of ion permeation enjoy an
important advantage over mean field models such as Gold-
man–Hodgkin–Katz theory (reviewed by Hille, 1992; re-
cently used by Dieckmann et al., 1999) or Poisson-Nernst-
Planck theory (for example, Chen et al., 1997; also
Kurnikova et al., 1999). In these models, a probability
distribution for ion concentration within the pore corre-
sponds to an average over states of 0, 1, 2, 3, . . . ions in the
pore.
Transitions between states of rate theory models are
exponentially distributed in time. When these transitions
describe ions within the pore, the exponential distributions
can be viewed as asymptotic approximations to diffusion
over energy barriers (Cooper et al., 1985, 1988). When the
transitions describe ions from the bulk solution entering the
pore, the exponential distribution corresponds to the as-
sumption that the ion entry rate into an empty channel does
not depend on the time elapsed since the channel last
became empty (McGill and Schumaker, 1996).
However, rate theory models are not entirely satisfactory
because they do not describe ion transport well when the
conditions for the asymptotic approximation for diffusion
over a barrier are not satisfied (Levitt, 1986; Dani and
Levitt, 1990). A way to overcome this difficulty is found in
the work of Levitt (1986), who showed how occupancy
restrictions can be incorporated into diffusion models.
McGill and Schumaker (1996) demonstrated that Levitt’s
model can be viewed as a diffusion within a state diagram
analogous to rate theory. The diagram is parameterized by a
single continuous reaction coordinate. Those authors further
showed how Levitt’s boundary conditions can be modified
so that diffusers entering the pore are exponentially distrib-
uted in time, corresponding to ions entering an empty chan-
nel at steady state.
The framework model we construct below is designed to
incorporate the molecular dynamics results of Pome`s and
Roux (1996, 1997, and manuscript in preparation), who
show how proton permeation may be dependent on both the
potential of mean force of an excess proton within the
permeation pore and the potential of mean force of water
Framework Model for H Conduction 13
Biophysical Journal 80(1) 12–30

reorientation in the empty pore (that is, without an excess
proton). In the first section, we construct the model as the
diffusion limit of a random walk. The purpose of this
construction is to obtain boundary conditions that restrict
pore occupancy to a single excess charge, or a single defect
in water orientation. The diffusion limit of a random walk is
very well known in both the physical (Chandrasekhar,
1943) and mathematical literature (e.g., Karlin and Taylor,
1981). A good introduction can be found in the first chapter
of Zauderer (1989). However, our discussion is self-con-
tained. The boundary conditions we obtain allow the result-
ing model to be solved analytically.
The model is then solved for the special case of thermo-
dynamic equilibrium, showing how the occupation proba-
bility of the pore is related to the corresponding expression
from statistical mechanics. We obtain the general solution
for the current through the framework model and discuss
several of its properties. The solution depends on potentials
of mean force and diffusion coefficients, which can be
obtained from molecular dynamics. We further obtain the
solution for the channel conductance, with the derivative of
current with respect to applied voltage evaluated at thermo-
dynamic equilibrium. Finally, we describe the method of
numerically simulating trajectories underlying the frame-
work model. Schumaker et al. (2000) incorporate into the
framework model the results of the molecular dynamics
simulations of Pome`s and Roux (manuscript in preparation).
They then make a detailed comparison with experiment.
CONSTRUCTION OF THE MODEL
Figure 1 A shows a state diagram for the dynamics of proton
permeation through gramicidin based on the results of
Pome`s and Roux (manuscript in preparation). The state
diagram consists of two segments. The top segment corre-
sponds to diffusion of an excess proton through the pore.
For simplicity, we will sometimes refer to these states as a
proton occupying the pore. The excess proton cannot be
uniquely identified. Pome`s and Roux used as their reaction
coordinate the axial component of the orientation moment
of the pore contents. Schumaker et al. (2000) show how this
may be rescaled to give the axial component of the dipole
moment of the pore contents computed with respect to an
origin at the center of the channel; we denote this quantity

H. To illustrate the meaning of H, consider a simple
example. Let z be the spatial coordinate co-axial with the
pore, extending over the interval L/2  z  L/2. If we
ignore pore waters and consider only an occupying proton,
then e0L/2  
H
 e0L/2.
The proton segment in Fig. 1 A is parameterized by H,
with values ranging over the interval [A
H, A
H]. The value

H
 A
H corresponds to a proton at the channel entrance
on the left (side I), and H  A
H corresponds to a proton at
the channel entrance on the right (side II). In general, we
expect A
H
e0L/2, since the polarization of the pore waters
in response to the excess charge will reduce the dipole
moment (Roux and Karplus, 1993). The cartoons at the
upper right and left hand corners of the diagram depict
states on the proton segment.
Figure 1 B shows the intrinsic potential of mean force,

H, computed for the proton reaction coordinate (Pome`s
and Roux, 1997, and manuscript in preparation). The sim-
ulated system included the channel, pore waters, and a few
waters clustered outside each channel entrance. The mem-
brane lipid and bulk aqueous solution were not simulated.
The word intrinsic (Levitt, 1986) refers to components of
the potential apart from that due to the experimentally
applied transmembrane potential. H has a shallow poten-
tial minimum near the center of the interval, corresponding
to the excess charge located near the center of the pore.
The bottom segment of Fig. 1 A is parameterized by the
dipole moment of the water molecules in the pore in the
absence of an excess proton, denoted d. For simplicity, we
will sometimes refer to these states as corresponding to an
empty pore or a defect occupying the pore, and we will refer
to the bottom segment as the defect segment. Values of d
range over the interval [A
d , A
d ]. The value d  A
d
corresponds to water dipole aligned with oxygens pointing
to the right and d  A
d corresponds to water dipoles
aligned with oxygens pointing to the left.
The molecular dynamics simulations suggest that diffu-
sion of this reaction coordinate is often associated with an
entrance-initiated defect in the hydrogen bond chain, with
water dipoles aligned on either side as suggested by the
cartoons below the defect segment. This terminology was
introduced by Phillips et al. (1999), and refers to a defect
that originates on the side of the channel opposite an exiting
proton. Those authors have suggested that defects may be
exit-initiated instead. Formally, the single proton conduc-
tion model that we develop here does not depend on this
choice.
Figure 1 C shows the intrinsic potential of mean force,

d, computed for the defect reaction coordinate. Note that
values on the abscissa increase from right to left. The
potential minima are reflected in the molecular dynamics
simulations by defects frequently found near one of the
channel entrances, as indicated by the cartoons at the lower
left and lower right of Fig. 1 A. These results suggest that a
proton may enter a channel that has d in a range of values
concentrated near one of the potential minima. This set of
possible transitions is suggested by the dashed lines be-
tween the segments in Fig. 1 A.
Figure 2 A again shows the state diagram corresponding
to the dynamics of proton permeation, with the set of
possible transitions between segments denoted by dashed
lines. These divide the defect segment into 3 regions; com-
pare with Fig. 1 C. The interior, corresponding to most of
the central barrier, occupies the interval [C
d , C
d ]. Sur-
rounding this are boundary region I to the left and boundary
14 Schumaker et al.
Biophysical Journal 80(1) 12–30

region II to the right. Region I is the interval of defect
reaction coordinate from which a proton can enter the
channel on side I, and region II is the corresponding interval
on side II.
The gramicidin dimer is physically symmetrical about the
center of the pore, and the intrinsic potentials calculated by
the molecular dynamics simulations are very nearly sym-
metrical about the midpoints s  0, s  {H, d}. We will
assume that these potentials are exactly symmetrical. Fur-
ther, there is an unknown energy difference between the
FIGURE 1 Hypothetical proton conduction mechanism. All energies are
in units of kBT for T  298 K. (A) State diagram for proton conduction
mechanism. The top segment is parameterized by the proton reaction coordi-
nate,H. Cartoons at upper left and upper right depict pore contents at the ends
of the segment. An excess proton may enter from side I, at the left, and pass
through the pore to exit on side II, at the right. Pore waters are depicted as
angles with oxygen at the vertex. Oxygens tend to align toward the proton. The
bottom segment is parameterized by the defect reaction coordinate, d. A
defect in the hydrogen bonding structure between waters must pass through the
channel so that waters are realigned to accept another proton from side I.
Dashed lines indicate that the transition from the proton-occupied state to the
defect state may occur for the defect in a range of locations. (B) Proton
potential of mean force, H (dots), and applied potential energy, H (solid).
The proton PMF shown is that calculated by Pome`s and Roux (1997).
Potentials are defined in the interval A
H

H
A
H. (C) Defect PMF, d
(dots), as calculated by Pome`s and Roux (1997, and manuscript in preparation)
and applied potential energy, d (solid). Energy minima at either end of the
central barrier correspond to a state with a defect near one end of the pore.
Intervals of reaction coordinate between dashed lines on either side of the
central barrier are the boundary regions. These are lumped together to form the
boundary states bI and bII, shown in Fig. 2 B.
FIGURE 2 State diagrams and random walks. (A) State diagram of
hypothetical proton conduction mechanism; compare with Fig. 1 A. A
H
and A
d are the extreme values of the reaction coordinate intervals;
compare with Fig. 1, B and C. The interior of the defect interval is bounded
by coordinates C
d . Boundary regions are the subintervals C
d

d

A
d . In the framework model they have effective electrical coordinates
B
d . The symbols f H and f A
d , f B
d, and f C
d denote the electrical width of
their respective subintervals. (B) State diagram of the framework model.
The boundary regions are lumped into boundary states bI and bII. (C)
Random walk used to construct the framework model. States H1, . . . , Hn
scale to the proton segment in the diffusion limit, n 3
. States d1, . . . ,
dn are ordered right to left and scale to the defect segment. States bI and bII
scale to the discrete boundary states. (D) Symmetrized random walk used
by the numerical simulations. The notation for transition probabilities is
similar to that shown in panel C.
Framework Model for H Conduction 15
Biophysical Journal 80(1) 12–30

proton and defect potentials of mean force shown in Fig. 1,
B and C. For these reasons, we will sometimes find it useful
to refer to the following relative potentials of mean force,

H

H

H

H
A
H , (1)

d

d

d

d
C
d , (2)
where A
H
 
H(A
H)  H(A
H), and C
d
 
d(C
d ) 

d(C
d ).
Applied Field
Let V(z) denote the component of the electrostatic potential
due to an applied transmembrane potential, VI, on side I. It
is independent of the charge distribution of the channel and
the pore waters (Roux, 1997). V(z) is assumed to drop
linearly over the length of the channel, corresponding to a
constant electric field E  VI/L in the positive z direction.
This form for the potential is appropriate for the simple
cylindrical geometry of gramicidin (Jordan et al., 1989;
Roux, 1999). Some feathering of the potential does occur
near the channel entrances in these calculations, but, for
simplicity, we will use the linear approximation.
When the applied electric field is constant, the resulting
contribution to the potential energy of the pore contents
depends only on its net charge and dipole moment (e.g.,
Jackson, 1975). This is seen by considering the electrical
potential energy, s, as a function of the charge density,

s(z), associated with species s and the electrical potential
V(z),

s


L/2
L/2

s
zVz dz, (3)
where the electrical potential is given by the linear drop
Vz VI/2 zE (4)
over the interval L/2 z L/2. For the cases of the
proton occupied and empty pores, the applied potential
energy can be written in the forms,

H
 e0VI/2HE, (5)

d

dE, (6)
where e0 is the elementary electronic charge, and 
s is the
dipole moment reaction coordinate. Electrical potential en-
ergies corresponding to VI  0.1 V are shown in Fig. 1, B
and C.
The potential energy drop around the cycle of the state
diagram in Fig. 2 A must be zero. To compute the drop,
consider the following energy differences:

H
A
H

H
A
H
 2A
HE, (7)

H
A
H

d
A
d
 0, (8)

d
A
d

d
A
d
 2A
dE, (9)

d
A
d

H
A
H
  e0VI , (10)
Eqs. 7 and 9 follow directly from Eqs. 5 and 6. To verify
Eq. 8, consider a proton leaving the channel at H  A
H on
side II. We assume that the orientation of the water dipoles
in this state and the state d  A
d are the same. Because
the electrical potential energy of the proton on side II is
zero, it then follows that H(A
H) and d(A
d ) are equal.
Similarly, Eq. 10 is obtained by considering that a proton
leaving the channel from side I carries with it a potential
energy e0VI. Adding Eqs. 7–10 and dividing by E gives
2A
H
 2A
d
 e0L. (11)
This equation relates the physical length of the pore with the
maximum values of the proton and defect reaction coordi-
nates. When the values A
H and A
d obtained from the
molecular dynamics simulations of Pome`s and Roux (manu-
script in preparation) are used, we obtain the result L  22.9
Å (Schumaker et al., 2000). This is slightly shorter than the
physical length of the pore. The discrepancy may be due to
the confinement of the excess proton in the molecular
dynamics simulations (Schumaker et al., 2000).
It is convenient to introduce dimensionless electrical dis-
tances, analogous to those encountered in rate theory (e.g.,
Hille, 1992).
f H  2 A
H
e0L, (12)
f A
d
  A
d
 B
d
/e0L, (13)
f B
d
  B
d
 C
d
/e0L, (14)
f d  f C
d
 2 C
d /e0L, (15)
where the notation f d is used in the appendices. These
electrical distances are proportional to widths of subinter-
vals of reaction coordinate as shown in Fig. 2 A. Eqs. 13 and
14 refer to the value B
d . The points d   B
d are shown
in Fig. 2 A and will be the effective electrical coordinates of
the boundary regions. Expressing Eq. 11 in terms of the
electrical distances, we have
2f A
d
 2f B
d
 f C
d
 f H  1. (16)
The total energy, W, of the pore contents is the sum of the
intrinsic and applied potentials,
Wssssss. (17)
Random walk limit to a diffusion process
In this section, we construct the framework model for pro-
ton diffusion through gramicidin, whose state diagram is
shown in Fig. 2 B. The random walk construction that we
use obtains the Smoluchowski equations (or their first inte-
16 Schumaker et al.
Biophysical Journal 80(1) 12–30

grals, the Nernst–Planck equations), which describe diffu-
sion of the reaction coordinates H and d in the proton-
occupied or proton-empty channels, respectively. Most
significantly, the construction also obtains the boundary
conditions that make possible a description of diffusion on
the state diagram, which is a simplified configuration space
for proton conduction through gramicidin. That is, by the
structure of the random walk, we describe either the diffu-
sion of a single excess charge through the channel or the
diffusion of the axial component of the dipole moment of
the pore waters in the absence of an excess charge.
The difference between the state diagram of Fig. 2 B and
that of Fig. 2 A lies in the description of the boundary
regions C
d

d
 A
d . The framework model lumps these
into discrete boundary states bI and bII in the lumped state
approximation. These states are constructed so that their
probabilities are equal to the integral of the Boltzmann
factor over the boundary regions under conditions of sym-
metrical equilibrium; see Eq. 64 below. This approximation
greatly simplifies the mathematical description of entrance
and exit. Instead of a continuum of possible transitions
between the boundary regions on the defect segment and the
endpoints of the proton segment, as suggested by Fig. 2 A,
we have a pair of transitions between the lumped states bI
and bII and the proton segment, shown in Fig. 2 B. The
lumped states surround the interior of the defect segment,
which contains the central barrier shown in Fig. 1 C. Trans-
port of the defect reaction coordinate d in the interior is
described diffusively, by a Nernst–Planck equation. This
lumped-state approximation gives an accurate description of
transport over the barrier (Schumaker et al., 2000; Mapes
and Schumaker, submitted) while leading to a model that is
analytically solvable.
The state diagram of the random walk is shown in Fig.
2 C. It is discrete in both space and time. State Hi, i  {1,
2, . . . , n}, denotes a proton at coordinate i
H
 A
H(2i/n 
1), and state di, i  {1, 2, . . . , n}, denotes a defect at
coordinate i
d
 C
d (2i/n  1). State si will refer to either Hi
or di. The two additional states bI and bII will be taken to the
boundary states of Fig. 2 B by the random walk construc-
tion. There are altogether 2n  2 states in the random walk.
We will define transition probabilities appropriately and
take the limit n 3
to obtain the framework model.
Nernst–Planck equation for channel interior
Let H  2A
H and d  2C
d be the lengths of the reaction
coordinate intervals over which transport will be described
by a diffusion process. The distance between the states si of
Fig. 2 C is s  s/n. Diffusion over the reaction coor-
dinate intervals is described by diffusion coefficients, s,
having units of (dipole moment)2/time; for simplicity, we
assume that these are constants, independent of s.
The probabilities, i
s, for a transition from state si to si1,
and i
s, for a transition from si to si1, are given by
i
s
 t

s

s

2 exp
1
2
Wsi
s
Wsi1
s


, (18)
i
s
 t

s

s

2 exp
1
2
Wsi
s
Wsi1
s


, (19)
where  (kBT)
1 and t is the time interval between steps
of this discrete-time random walk. We scale t with n so
that the leading order of n in the expressions for i
s and i
s
is n0; transition probabilities then remain positive and finite
in the limit n 3
. Let
t /n2, (20)
where  is independent of n.  must be chosen so that the
probability of leaving a state at each time step is no greater
than 1. This choice is made explicitly by the algorithm
described in Numerical Solution below.
The transition probabilities i
s and i
s lead to the Boltz-
mann distribution at equilibrium. To see this, let Qi
s be the
probability that state si of the random walk is occupied. At
equilibrium, the system is in detailed balance with zero net
flux between any two states. This means
Qi
s
i
s
 Qi1
s
i1
s . (21)
Inserting the definitions for the transition probabilities gives
the result expected from the Boltzmann distribution,
Qi
s
Qi1
s
 exp Wsi1
s
Wsi
s
. (22)
The transition probabilities may be expanded in the small
parameter 1/n to give a useful expression in preparation for
taking the limit n 3
. The expansion gives
i
s
 t

s

s

2 1
1
2

sWsi
s

 
s

2
i
s
 n3

, (23)
i
s
 t

s

s

2 1
1
2

sWsi
s

 
s

2
i
s
 n3

, (24)
where
i
s
Wsi
s
/4 2Wsi
s

2/8, (25)
and primes denote derivatives with respect to the argument.
We assume Ws is continuous. The transition probabilities
used by McGill and Schumaker (1996) are obtained from
Eqs. 23 and 24 by truncating after the first-order terms.
We now construct a limit of the random walk that leads
to the Nernst–Planck equation at steady state. Consider the
Framework Model for H Conduction 17
Biophysical Journal 80(1) 12–30

states si, 2  i  n  1. At steady state, the total probability
flowing into si at each time step equals the total probability
flowing out. Equating these flows, we have the expression
for probability balance,
Qi
s
i
s
 i
s
 Qi1
s
i1
s
 Qi1
s
i1
s . (26)
Each state, si, represents a segment of the s interval of
length s  s/n. While taking the limit n3
, we wish
to consider the density, Ps, related to the state probabilities
by
Qi
s
 Pi
s

s, where Pi
s
 Psi
s
. (27)
Substitute Eqs. 23, 24, and 27 into 26, simplify, and use the
following finite difference expressions for first and second
derivatives:
WssPss lim
n3

(Wi1
s
 Pi1
s
Wi1
s
 Pi1
s )/2s,
(28)
Pss lim
n3

(Pi1
s
 2Pi
s
 Pi1
s )/s2, (29)
where Wi
s
  Ws(i
s) and i
s
3 
s as n3
. Note from Eq.
25 that i
s tends to a continuous function of s in this limit.
We obtain the Smoluchowski equation,
0
d2Ps
ds2

d
ds
WsPs. (30)
This is a diffusion equation with the second term corre-
sponding to a systematic force on the diffuser whose mag-
nitude is proportional to the gradient of the potential energy
W. Integrating once, we obtain the Nernst–Planck equation,
Js s

dPs
ds
 WsPs

, (31)
where Js is the flux of species s. The first term on the
right-hand side is Fick’s law, and the second term corre-
sponds to the systematic force.
Note that Js is positive when it is in the direction of
increasing s. That means that the flux of protons is positive
when directed from left to right in Fig. 2 A, whereas the flux
of defects is positive when directed from right to left. In
both cases, positive flux corresponds to progress in the
clockwise direction around the diagram.
Entrance and exit transition probabilities
In this subsection, we obtain expressions for the entrance
and exit transition probabilities that connect the following
states in Fig. 2 C: H1, bI, and dn on side I of the pore and Hn,
bII, and d1 on side II. In particular, these probabilities are
scaled with n to obtain, in the limit n3
, the state diagram
of Fig. 2 B. In this figure, the interior of the defect segment
supports a probability density and the endpoints, bI and bII,
have positive probability. The limit on n is taken in the
Boundary conditions subsection, following this one.
We begin the formal development by considering proton
transitions into and out of the channel for the random walk
state diagram shown in Fig. 2 C. At equilibrium, there is no
net flow of ions into the channel at either entrance. On side
I, we form the detailed balance relationships between H1
and bI and between bI and dn and use these to eliminate QI
b;
the probability that boundary state bI is occupied. Similarly,
we eliminate QII
b on side II. This leads to
Q1
H
Qn
d 
ICI
I

I
I
, (32)
Qn
H
Q1
d 
IICII
II

II
II
, (33)
where the CR are the excess proton concentrations on side R.
We use concentrations instead of activities because diffu-
sion coefficients are defined in terms of concentration gra-
dients (Robinson and Stokes, 1965; see also McGill and
Schumaker, 1996).
We require that the detailed balance equations remain
satisfied in the limit n 3
. In this limit, Q1
H and Qn
H
converge to values proportional to the probability density
PH on the proton segment at the endpoints A
H and A
H,
respectively. Similarly, Q1
d and Qn
d converge to values pro-
portional to Pd(C
d ) and Pd(C
d ), respectively. The equi-
librium densities on the proton and defect segments are
obtained by solving Eq. 31 with Js  0, giving the Boltz-
mann distribution,
Pss KseWs(s), (34)
where the Ks are constants. Take the limit of Eqs. 32 and 33
as n 3
, using Eqs. 27 and 34. We obtain
lim
n3

ICI
I

I
I

KH
Kd
expWdC
d
WHA
H
, (35)
lim
n3

IICII
II

II
II

KH
Kd
expWdC
d
WHA
H
. (36)
From the left-hand side of these equations, we see that the
constant KH/Kd must be proportional to CI in the first of
these equations and to CII in the second. This may be
understood by considering the Nernst equation,
CII  CIeI, (37)
where I  eVI. We will set
KH
Kd

aˆCI
aˆC0
eI 
aˆCII
aˆC0
, (38)
where the dimensionless quantity, aˆ, and the concentration,
C0, have been introduced. We assume these do not depend
on CI, CII, or I. Below, aˆ is determined by the requirement
18 Schumaker et al.
Biophysical Journal 80(1) 12–30

that occupation probabilities of the framework model agree
with statistical mechanics in the case of a symmetrical
equilibrium (VI  0). We discuss C0 further in Equilibrium
Probability for Proton Occupation.
Insert the second expression of Eq. 38 into the right-hand
side of Eq. 35 and the third expression of Eq. 38 into the
right-hand side of Eq. 36. The resulting detailed balance
relationships are satisfied by the following decomposition
lim
n3

nICI
I

1
aˆ
CI
C0
expWdB
d
WHA
H
I,
(39)
lim
n3


I
n I
 aˆ expWdC
d
WdB
d
, (40)
lim
n3

nIICII
II

1
aˆ
CII
C0
expWdB
d
WHA
H
, (41)
lim
n3


II
n II
 aˆ expWdC
d
WdB
d
. (42)
In making this decomposition, the boundary points bI and
bII are formally assigned defect reaction coordinates B
d .
The additional term I in the exponent on the right-hand
side of Eq. 39 represents the electrostatic energy of an ion
entering on side I. The distribution of factors of n on the
left-hand side of Eqs. 39–42 reflects the fact that the rates
of transitions from the proton and defect segments into the
boundary states bI and bII must scale with one power of n
higher than the rates of transitions from the boundary states
back to the segments. This is because the boundary state
probabilities QI
b and QII
b remain positive in the limit n3

while the states s1 and sn, s  {H, d}, scale to the endpoints
of probability densities.
We next introduce the new quantities
a aˆ expC
d
B
d
, (43)
 C
d
A
H
 kBT ln C0/C• , (44)
where B
d
 
d(B
d )  d(B
d ) and C• is the unit con-
centration, e.g., C•  1 M (the argument of the logarithm
must be dimensionless). In the expression for , the term
C
d
 A
H depends on the absolute energy difference be-
tween the proton and defect potentials of mean force. This
energy difference was not determined by the molecular
dynamics. As a consequence,  will be treated as an adjust-
able parameter in our analysis of the Eisenman et al. (1980)
conductance data using the single proton model (Schumaker
et al., 2000).
We now simplify the exponents of Eqs. 39–42 by replac-
ing aˆ with a, decomposing W according to Eq. 17, express-
ing H in terms of d using Eqs. 8 and 10, using the
definitions for electrical distances, Eqs. 13 and 14, and
finally introducing the definition of . The result is
lim
n3

nICI
I

1
a
CI
C•
expf A
d
I  , (45)
lim
n3


I
n I
 a expf B
d
I, (46)
lim
n3

nIICII
II

1
a
CII
C•
expf A
d
I  , (47)
lim
n3


II
n II
 a exp f B
d
I. (48)
The following definitions are consistent with these con-
straints:
I  II  tta1a1C•
1e, (49)
I  tta1nexp f A
d
I , (50)
II  tta1nexpf A
d
I , (51)

I  ttd1n2aexpf B
d
I , (52)

II  ttd1n2aexp f B
d
I , (53)
I  II  tt d1n, (54)
where the access time, ta, is the second free parameter
optimized to fit the Eisenman et al. data set. We also
introduce the proton and defect diffusion times
ts  s2/s, (55)
where s  {H, d}. Eqs. 49–54 complete our description of
the discrete-time, discrete-state random walk whose state
diagram is shown in Fig. 2 C. The time step is t. The
transition probabilities RCR, R,
R, and R, R  {I, II},
are all dimensionless. The rates of transitions into the proton
segment, equal to RCR/t, are independent of n. The
scaling of the R with n is then determined by the con-
straints given in Eqs. 45 and 47.
Our scaling of the entrance rate as independent of n
means that proton entrances are exponentially distributed in
time when the channel is accessible to protons. For exam-
ple, when the system is near boundary state bI of Fig. 2 B,
one waits an exponentially distributed length of the accu-
mulated time spent in state bI before a proton enters the
channel. This is very similar to rate theories, where transi-
tions between the discrete states of the system are exponen-
tially distributed in time and, in particular, one waits an
exponentially distributed time in the empty state before an
ion enters the channel. The exponential distribution may be
considered as arising from the assumption that the ion entry
rate into an empty channel does not depend on the time
elapsed since the channel last became empty (McGill and
Schumaker, 1996). This will often be a good approximation
for the entry of distinct ions into the channel, although it
Framework Model for H Conduction 19
Biophysical Journal 80(1) 12–30

does not capture the fractal character of the trajectories of
individual ions.
The electrical distance, fA
d , appearing in Eqs. 45 and 47 is
assigned to the R. This means that the entrance rate is
independent of I. Physically, this choice assumes that
entrance is limited by diffusion up to the channel mouth,
and not by a local barrier at the channel entrance. Diffusion
up to the channel mouth should not be very sensitive to an
applied electric field, which would not penetrate far into an
electrically conducting solution. In contrast, an applied field
might well speed transitions over a local barrier near the
membrane surface. The sublinear IV curves observed in the
proton conduction experiments of Eisenman et al. (1980),
Decker and Levitt (1988), and Akeson and Deamer (1991)
at low proton concentrations, where conductance is likely to
be limited by entrance, are consistent with the idea that the
entrance rate is not strongly dependent on I. Even in the
absence of a local entrance barrier, however, our assignment
of fA
d neglects interfacial polarization (Andersen, 1983),
which becomes significant at low ionic strength and large
applied potentials.
The transition probabilities,
R, are scaled with n in the
same way as i
s and i
s (see Eqs. 23 and 24, note s  1/n).
These latter transition probabilities lead to a diffusion pro-
cess in the limit n 3
. The scaling of the
R is thus
consistent with modeling motion of defects from the region
of the central barrier to the boundary regions as diffusive.
The scaling of the R with n is then determined by the
constraints given in Eqs. 46 and 48. The assignment of the
electrical distances, fB
d , to the transition probabilities,
R,
instead of the R is arbitrary. But we do not expect these
rapid transitions to be rate limiting, so this choice may not
be important.
Boundary conditions
In this subsection, we start from expressions for a steady-
state flux through the entrances I and II of the pore as
described by the random walk in Fig. 2 C. Inserting the
definitions for the transition probabilities that we have pre-
viously obtained, and taking the limit n 3
, we obtain
boundary conditions for the framework model of Fig. 2 B.
The calculations are straightforward, but somewhat cum-
bersome. Use of the dimensionless variables described in
Appendix A shortens the intermediate expressions.
On side I, balance probability flux in and out of the state
H1 of Fig. 2 C. We have
Q1
H
1
H
 I QI
b
ICI  Q2
H
2
H , (56)
Substitute in expressions for transition probabilities, Eqs.
23, 24, 49, and 50, as well as the relationship between state
probability and density, Eq. 27. Divide by a common factor,
rearrange terms, take the limit n 3
, and, finally, use the
definition of JH, Eq. 31, to obtain (see Appendix A)
PHA
H

H exp f A
d
IJHta  QI
ba1ecI ,
(57)
where cI  CI/C• is the dimensionless concentration.
To complete the analysis of side I, balance flux in and out
of state dn,
Qn
d
n
d

I QI
b
I  Qn1
d
n1
d . (58)
Make the appropriate substitutions, rearrange terms and
take the limit n 3
to obtain
QI
b
 PdC
d

da expf B
d
I. (59)
Similarly, on side II, flux in and out of state Hn is balanced
to obtain
PHA
H

H expf A
d
I JHta  Q II
ba1ecII , (60)
and balancing flux in and out of the state d1 gives
QII
b
 PdC
d

daexp f B
d
I. (61)
We can further eliminate QI
b and QII
b from Eqs. 57 and
59–61 to yield
PHA
H)Hexp f A
d
I
JHta  Pd C
d

dexpf B
d
I  cI , (62)
PHA
H

Hexpf A
d
I
 JHta  PdC
d

dexp f B
d
I  cII , (63)
Note that these expressions are dimensionless. The density,
Ps, has units of inverse reaction coordinate, and the flux, Js,
has units of inverse time.
Expression for a
The dimensionless quantity aˆ was originally introduced in
Eq. 38 and then expressed in terms of a in Eq. 43. To obtain
an expression for a, we set the probability of the boundary
states bI and bII equal to the integral of the Boltzmann factor
over the respective boundary regions under the condition of
symmetrical equilibrium: CI  CII and I  0. For example
QI
b


C
d
A
d
Pdd dd. (64)
At equilibrium, Pd(d) is given by Eq. 34. Substitute in for
QI
b using Eq. 59, with I  0, and then on both sides for
20 Schumaker et al.
Biophysical Journal 80(1) 12–30

Pd(d) using Eq. 34. Solving for a, we obtain
a d1

C
d
A
d
exp dd dd
 
d

1

A
d
C
d
exp dd dd. (65)
The value of a is given as an integral over the relative
energy, (d)  (d)  (C
d ), and can be evaluated
after a choice is made for C
d . Eq. 65 completes the con-
struction of the single proton model using the lumped state
approximation. We next show that the model is consistent
with statistical mechanics under the condition of symmet-
rical equilibrium.
EQUILIBRIUM PROBABILITY FOR
PROTON OCCUPATION
Insight may be gained into the physical assumptions that lie
behind the framework model by comparing the expression
for the equilibrium probability of proton occupation of the
pore with the analogous expression from statistical mechanics.
Statistical mechanics
For the purpose of constructing a partition function, the total
energy of the empty pore includes the energy of a reference
proton in solution far away from the channel. Assuming that
the reference state is on side I, and approximating activities
by concentrations, the energy of the reference proton may
be written in the form
Wref
H
 kBT ln CI/CH I , (66)
where CH is a constant with units of concentration. Accord-
ing to statistical mechanics, the probability of proton occu-
pation is then equal to
QSM
H


A
H
A
H exp WHH dH

A
H
A
H exp WHH dH
 CH/CIexp I
A
H
A
H exp W dd dd.
(67)
Because the intervals [A
H, A
H] and [A
d , A
d ] only in-
clude the pore interior, this expression neglects all interac-
tions between the channel and ions in the aqueous solution
outside.
Framework model
The equilibrium probability for proton occupation in the
framework model is
QFW
H


A
H
A
H
PHH dH. (68)
We also have the normalization condition
1

A
H
A
H
PHHdH


C
d
C
d
Pdddd  QI
b
 QII
b . (69)
We construct an expression comparable to that from
statistical mechanics by dividing the right-hand side of Eq.
68 by the right-hand side of Eq. 69. Form the ratio and
simplify (see Appendix B) to obtain
QFW
H


A
H
A
H exp WHH dH

A
H
A
H exp WHH dH 

H

d
C0
CI
exp I
,
(70)
where


C
d
C
d
exp Wdd dd
 exp dB
d


C
d
A
d
exp dd dd
 exp dB
d


A
d
C
d
exp dd dd.
(71)
To make the correspondence between QFW
H and QSM
H ,
equate the coefficient of the integral in the second term of
the denominator of Eq. 67 with the coefficient of  in Eq.
70. This leads to
CH  H/dC0 . (72)
With this identification, the equilibrium expressions for
QSM
H and QFW
H are identical when I  0. Under the influ-
ence of an applied potential, the framework model makes
the approximation

A
d
A
d
exp Wdddd . (73)
Framework Model for H Conduction 21
Biophysical Journal 80(1) 12–30

The expression for , Eq. 71, shows that the boundary
regions of the framework model act electrically like points.
C0 is a free parameter in the framework model theory that
we have developed thus far, and is combined with the
unknown free energy difference, C
d
 A
H, between the
defect and proton potentials of mean force in the expression
for  given by Eq. 44. Our analysis of the data of Eisenman
et al. (Schumaker et al., 2000) only requires a value for .
However, an expression for C0 can be obtained with the
help of the statistical mechanical theory of selective ion
channels developed by Roux (1999). This gives C0 as the
inverse of an effective pore volume accessible to the excess
proton (Schumaker, manuscript in preparation). This result
may be useful for interpreting temperature dependence stud-
ies of proton conduction.
GENERAL SOLUTION OF FRAMEWORK MODEL
In Appendix C, we find the general solution for the frame-
work model proton current, JH. The method of solution
depends on solving a linear set of equations, and, in this
respect, is similar to that of rate theory. See, for example,
La¨uger (1973) or Hille and Schwarz (1978).
A linear system of eight equations and eight unknowns is
obtained. Two of the equations arise by integrating the
Nernst–Planck equations, Eq. 31, for s  {H, d}. Four
equations come from the boundary conditions, namely Eqs.
57 and 59–61. The seventh equation is the normalization
condition, Eq. 69, and the eighth equation is simply the
equality of the proton and defect fluxes, JH  Jd, which
holds at steady state. The eight corresponding unknowns are
the densities at either end of the proton and defect diffusion
intervals, namely PH(A
H), PH(A
H), Pd(C
d ), and Pd(C
d ),
the probabilities QI
b and QII
b of the defect boundary regions,
and the fluxes JH and Jd. Three variables are quickly elim-
inated to give a system of five equations in Appendix C,
which is then solved.
The solution for the current has the form
JHcI , cII , I
cIeI  cII
tHFH  tdFd  taFa
, (74)
where the dimensionless variables cI and cII give the nu-
merical value of the bulk concentrations on either side of the
channel in standard units (e.g., moles), and I  eVI/kBT
is proportional to the applied electrical potential, VI. t
H and
td are the characteristic times for proton and defect diffusion
defined by Eq. 55, and the access time, ta, is a characteristic
time for entrance and exit, processes that were not repre-
sented in the molecular dynamics simulations. The coeffi-
cients of the characteristic times have the form
FH  F0
H
 F1
HcI  F2
HcII , (75)
Fd  F1
dcI  F2
dcII  F3
dcIcII , (76)
Fa  F0
a
 F1
acI  F2
acII , (77)
where the Fi
s are coefficients given in Appendix C.
Note that the coefficient Fd in the denominator of JH
includes a term proportional to the product cIcII. As cI and
cII increase, this term, in principle, gives an eventually
decreasing flux, JH. Single-ion theories that describe the
empty channel by a single indivisible state give rise to
expressions for JH whose forms are similar to Eq. 74, but
without such quadratic terms in the denominator. This is
true for single-ion rate theories (for example, La¨uger, 1973)
and single-ion diffusion theories (Levitt, 1986; McGill and
Schumaker, 1996).
The origin of the quadratic term is reminiscent of the
“clogging” mechanism described by Schumaker and Mac-
Kinnon (1990). When both cI and cII are sufficiently high,
the system becomes trapped in the proton segment of the
state diagram of Fig. 2 B. For example, an ion that exits on
side II will leave the system in state bII, but an ion entry
from side II is likely to occur before the defect reaction
coordinate can cross the central barrier to state bI. This
effect reduces the rate of cycling around the diagram and so
reduces the flux JH at very high cIcII.
PROPERTIES OF THE GENERAL SOLUTION
Current symmetry
The gramicidin dimer is symmetrical about its center, z  0,
a property reflected in the symmetry of the potentials of
mean force calculated by Pome`s and Roux (manuscript in
preparation): s(s)  s(s). As a consequence of this
and the symmetry properties of the applied potential 
about the channel center, the framework model current J,
Eq. 74, should have the symmetry
Jc2 , c1 , Jc1 , c2 ,  (78)
for all concentrations c1 and c2, as well as electrical poten-
tial energies. Further, the symmetry should hold while the
characteristic times, ts, are varied independently of each
other. Write the coefficients F i
s in the denominator of J,
given by Eqs. 75–77, as functions of the applied potential:
F i
s
 F i
s (I). Then the current symmetry holds if and only
if the coefficients have the properties,
F 0
s
 eF 0
s
, (79)
F 2
s
 eF 1
s
, (80)
F 3
s
 eF 3
s
, (81)
In Appendix D, we demonstrate that these coefficient sym-
metries hold.
22 Schumaker et al.
Biophysical Journal 80(1) 12–30

Saturating current at high voltage
At very high applied potentials, it is reasonable to expect
that the proton current through gramicidin should saturate at
a finite value. This is because one would not expect the
applied potential to reach far into the aqueous solutions on
either side of the membrane, because these solutions are
electrical conductors. This reasoning is consistent with IV
measurements through gramicidin by Akeson and Deamer
(1991). However, it ignores the effects of interfacial polar-
ization, which were measured by Andersen (1983) for the
conduction of several different cations through gramicidin
at low ionic strength.
The framework model is constructed to give the saturat-
ing behavior; the entrance rate constants, RCR, are inde-
pendent of the applied potential I, see Eq. 49. In Appendix
E, the limiting current is calculated to be
lim
I3

J cItaa1e. (82)
This result is easily understood in view of the expression for
RCR given by Eq. 49. The saturating current is simply the
entrance rate: the transition probability divided by the time
step t.
Equilibrium conductance
The conductance at Nernst equilibrium is defined by
GEq 
dI
dVI

Eq
 e0
2
dJH
dI
, (83)
where I  e0J
H is the current through the channel. Appendix
F shows how this quantity is computed from Eq. 74 and that
the denominator can be written in a surprisingly compact
form. The result simplifies even further under the condition
of a symmetrical equilibrium, with cI  cII  c and I 
0. Denoting by G0 the conductance evaluated at symmetri-
cal equilibrium, we find
G0 
e0
2
kBT
c
tHh0
H
 2tae  tdh0
dcg1
d
 g0
Hec
, (84)
where
h0
H
 
H

1

A
H
A
H
expHd, (85)
h0
d
 
d

1

C
d
C
d
expdd, (86)
g0
H
 
H

1

A
H
A
H
exp Hd, (87)
g1
d
 
d

1

A
d
A
d
exp d()) d. (88)
The first two terms in the expression for , Eq. 44, give
the free energy difference between the proton and defect
potentials of mean force. The denominator of Eq. 84 con-
tains terms proportional to both e and e. Consequently,
the conductance decreases if the free energy difference
between the profiles is either sufficiently large and positive
or sufficiently large and negative.
The conductance also decreases if the s have shapes
with high maxima (the h0
s are large) or deep minima (the g0
s
are large). Conductance is maximum when the s are
relatively flat and the magnitude of the energy difference, ,
is not too large. Finally, notice that the denominator of G0
is quadratic in c, reflecting the quadratic dependence of the
denominator of J, Eq. 74, on concentration.
NUMERICAL SOLUTION
Numerical simulations of random walks are used by Schu-
maker et al. (2000) to qualitatively demonstrate the nature
of the trajectories underlying the single proton conduction
model. They may also represent an interesting alternative
method for generating Brownian dynamics simulations (for
example, Jakobsson and Chiu, 1987). In this section, we
describe how to construct simulations corresponding to the
analytical theory.
The random walk shown in Fig. 2 C is convenient for
construction of the framework model, but less suitable for
numerical simulation. Instead, the numerical simulation
takes place on the symmetrized random walk shown in Fig.
2 D. To account for the additional state, we replace n3 n 
1 in the formula for the time step, Eq. 20, and in the
formulas for the entrance and exit transition probabilities,
Eqs. 49–54. These formulas then define the simulated ran-
dom walk.
According to Eq. 27, each random walk gridpoint corre-
sponds to a reaction coordinate subinterval of width s.
This is true even for the symmetrized endpoints i  0 and
i  n corresponding to coordinates d  C
d . The sym-
metrized random walk should then most appropriately be
identified with diffusion on the interval (1  1/n)C
d . For
this reason, we compare analytical results using the value
for a given by Eq. 65 with numerical results using the
modified value
asim  d)1

(11/n)C
d
A
d
exp(d(d)) dd. (89)
In the simulations, half of the probabilities associated with
the endpoints i  0 and i  n are added to the boundary
state probabilities QI
b and QII
b , respectively. With our defi-
Framework Model for H Conduction 23
Biophysical Journal 80(1) 12–30

nition of the symmetrized random walk, use of asim and this
treatment of probability is necessary for the good agreement
between the simulated and exact probability densities de-
scribed below. However, this method of simulating the
boundary conditions has an ad hoc character and can un-
doubtedly be improved.
The simulated random walk is discrete in both time and
space, conforming to the theoretical discussion. Because of
the large difference between the proton and defect diffusion
coefficients obtained from molecular dynamics (Schumaker
et al., 2000), transition rates on the proton segment of the
state diagram are much higher than those on the defect
segment. To retain numerical efficiency, different values of
 are used to describe proton and defect diffusion. For the
same reason, dwell times in the states bI and bII are calcu-
lated from geometric distributions directly, without taking
time steps, similar to the method used to calculate exit from
the empty state of the single ion model (McGill and Schu-
maker, 1996).
Proton and defect probability densities are shown in Fig.
3, A and B. These runs use the potentials of mean force and
diffusion coefficients obtained from the molecular dynam-
ics simulations of Pome`s and Roux (manuscript in prepa-
ration), and shown in Fig. 1, B and C. Densities shown in
Fig. 3 represent the average of four simulations of about
2.1  109 time steps each. Each of the four simulations was
given a different seed to initialize the random number
generator (ran2, described by Press et al., 1992) but all other
parameter values were kept the same. In making these
calculations, the intrinsic potentials are interpolated linearly
between grid points where potentials of mean force were
calculated by molecular dynamics. The piecewise smooth
character of the probability densities is due to the piece-wise
linearity of the interpolated potentials.
The simulated and exact proton densities are very similar
in both shape and amplitude. The integrated probability of
the simulated proton density is only about 0.5% greater than
that of the exact solution for the proton density. This excel-
lent agreement is in part fortuitous, because only 896 pro-
tons entered and left the pore in the simulations that pro-
duced Fig. 3. If we estimate the uncertainty in this number
as 8961/2, the corresponding variation in the probability for
proton occupation is about 3.3%. In a second simulation,
made using a different parameter set, the difference between
the simulated and exact proton occupation probability was
slightly larger than the estimate based on the square root of
the number of protons that enter the pore (results not
shown). For the simulation used to make Fig. 3, the defect
density and boundary state probabilities QI
b and QII
b are also
in excellent agreement with the exact solution.
Fig. 4 shows the current through the framework model as
a function of applied voltage. The exact solution of the
framework model is compared with simulations at 50, 100,
and 200 mV applied potential. Agreement is very good. We
have also used the numerical simulation method to calculate
mean first passage times across the defect potential barrier
as a function of applied potential, following the protocol
used to construct Fig. 5 in Schumaker et al. (2000). Mean
first passage times obtained by the simulation method are in
excellent agreement with those results (unpublished).
DISCUSSION
We have constructed a framework model for single proton
conductance through gramicidin, a model designed to in-
corporate the results of molecular dynamics simulations
(Pome`s and Roux, 1996, 1997, and manuscript in prepara-
tion) and used it to predict conductance properties. Such
models are currently necessary to bridge the gap in time
FIGURE 3 Exact and simulated framework model probability densities.
The analytical solution of the framework model is given by the solid
curves, and results of numerical simulations are given by the dots. Proton
and defect trajectories were simulated on the symmetrized random walk
shown in Fig. 2 D with n  57. Densities were obtained by averaging
together the results of four runs of about 2.1  109 time steps each. CI 
CII  0.01 M and VI  100 mV. Parameter values are A
H
 3.35, A
d

8.1, B
d
 6.8, and C
d
 5.7 with   0.4348 and scaling d  d (see
Schumaker et al., 2000). In addition,  3.60 and ta  23.5 ns. (A) Density
on the proton segment for A
H

H
A
H. (B) Density on the defect
segment for C
d

d
C
d . The analytical solutions for the probabilities
of the boundary regions on either side of the defect interval are QI
b
 0.440
and QII
b
 0.188. Numerically estimated values of these probabilities are
QI
b
 0.438 and QII
b
 0.184.
24 Schumaker et al.
Biophysical Journal 80(1) 12–30

scales between molecular dynamics simulations and elec-
trophysiology. The framework model consists of a coupled
pair of Nernst–Planck equations describing the diffusion of
the axial component of the pore dipole moment in the
presence or in the absence of an excess proton in the chain
of water molecules occupying the gramicidin pore. In the
absence of an excess proton, the dipole moment parameter-
izes the reorientation of the pore waters from a configura-
tion in which the axial component of their dipole moments
point in the z direction to a configuration in which the
axial component points in the z direction. The simulations
suggest that water reorientation is mediated by a defect in
the hydrogen bond structure that passes through the water
chain. Therefore, we also describe the empty or neutral
channel as one in which a defect occupies the pore.
The Nernst–Planck equations and their boundary condi-
tions are obtained as the limit of a sequence of random
walks. By the random walk construction, the boundary
conditions restrict the model to describing a single diffusing
proton or a single diffusing defect. In other words, we
construct a diffusion process on a simplified configuration
space for proton conduction through gramicidin, portrayed
by the state diagram shown in Fig. 2 B. This is much
different from Goldman–Hodgkin–Katz theory (Hille,
1992) or Poisson Nernst–Planck theories (for example,
Chen, 1997), which are mean field theories. In these latter
cases, a solution corresponds to a superposition of many
different occupancy states of ions in the pore.
From the simulations, we obtain the potential of mean
force for proton and defect diffusion through the pore as
well as estimates of the proton and defect diffusion coeffi-
cients (Schumaker et al., 2000). The proton potential of
mean force, shown in Fig. 1 B, has the form of a shallow
potential well. Transport across such a profile is not de-
scribed well by rate theories, which are based on asymptotic
approximations to diffusion over energy barriers (Cooper et
al., 1985, 1988). The defect potential of mean force has the
form of a high central barrier separating potential minima.
The potentials and diffusion coefficients, obtained from
molecular dynamics, completely determine the description
of transport in the pore interior. By incorporating this in-
formation from the simulation, the framework model can
describe this portion of the permeation process without
adjustable parameters.
However, the simulations of Pome`s and Roux do not
address the processes of ion entrance and exit. Here, it is
necessary for the framework model to introduce some de-
scription. We take a very simple and phenomenological
approach. Proton entrances on a given side are assumed to
be exponentially distributed in time when the channel is
accessible to occupation. The entrance process is then de-
scribed by a single adjustable parameter, the mean time
before a proton enters an empty channel. This is the same
description of ion entrance as is implicit in rate theories. The
exit process must then satisfy detailed balance with entrance
at equilibrium. It also involves a single adjustable parame-
ter, which would be known if the absolute free energy
difference between the proton and defect potentials of mean
force were known. Finally, there is a third adjustable pa-
rameter, which is also concerned with the ion entrance and
exit process. This is due to the fact that the state of the
waters in the empty channel is not precisely known just as
an excess proton leaves or enters; instead, the dipole mo-
ment of the empty channel may be in one of the boundary
regions of Fig. 1 C. The width of these regions is the third
parameter. In practice, a reasonable value for this parameter
can be chosen, encompassing the potential minima on either
side of the central barrier (Schumaker, et al. 2000). Only the
two parameters controlling the entrance and exit rates need
be optimized to make a comparison with experiment.
The conduction mechanism of our framework model is
easily grasped by means of a state diagram, which is anal-
ogous to those that describe rate theory permeation models,
see Fig. 2 B. As a proton passes through the channel, say
from side I to side II, the proton reaction coordinate follows
its progress on the upper segment of the diagram from left
to right. The proton exits the channel, leaving water mole-
cules partially aligned. This corresponds to the lumped state
bII. The dipole moment of the water chain must then turn to
become receptive to another proton entering on side I. The
molecular dynamics simulations (Pome`s and Roux, 1997,
and manuscript in preparation) suggest that turning is me-
diated by a defect in the hydrogen bond structure that passes
through the water chain. This defect is followed by the
reaction coordinate that parameterizes the lower segment of
the state diagram, from right to left. When the reaction
coordinate is in the lumped state bI, the system is ready to
accept another free proton from side I.
The most problematic feature of this model may be the
treatment of the boundary conditions. An important com-
FIGURE 4 Exact and simulated framework model proton current. The
figure shows current I as a function of applied voltage under the same
conditions as in Fig. 3. The analytical solution of the framework model, Eq.
74, is shown by the solid curve. Error bars show twice the standard
deviation of the mean for the four runs described in the legend of Fig. 3.
Framework Model for H Conduction 25
Biophysical Journal 80(1) 12–30

ponent of the configuration space is missing, which was
pointed out following the statistical mechanical expression
for the probability of proton occupation, Eq. 67. Therefore,
there is no description of the interaction of the channel with
excess protons outside, but close to, the entrance. For ex-
ample, it is plausible that water dipoles are partially oriented
by nearby protons before they enter the channel. In contrast,
our model implicitly assumes that the orientation of channel
waters is independent of excess protons on the outside, and
that proton entrance only becomes possible when the water
dipoles happen to be favorably aligned. To test this assump-
tion and possibly construct a more realistic model, it would
be very useful to have molecular dynamics simulations of
the ion entrance process.
APPENDIX A: DIMENSIONLESS VARIABLES AND
BOUNDARY CONDITIONS
Results in the body of this paper are given in dimensional form, but many
calculations become less cumbersome when dimensionless variables are
used. Reaction coordinate intervals over which proton or defect diffusion
is described are conveniently rescaled to [0, 1]. Let

H
 1H/A
H
/2, (90)

d
 1d/C
d
/2. (91)
Then the random walk of Fig. 2 C takes place on the grid points i
s
 i/n,
{i  1, 2, . . . , n}.
Dimensionless energies are defined: s(s)  s(s), s(s) 

s(s),
s(s)  s(s), and ws(s)  Ws(s). For the energy due to
the applied potential, write
I  I. The dimensionless probability
density is ps(s)  Ps(s)s. Electrical potential energies, corresponding to
Eqs. 5 and 6, can then be written:

H

H

I/2 fH2H  1
I/2, (92)

d

d
f C
d
2d  1
I/2. (93)
For the derivation of the boundary conditions, it is helpful to use the
compact notation, wi
s
 ws(i
s) and wi
s
  ws(i
s), where the prime denotes
a derivative with respect to the argument. The transition probabilities can
then be written as
i
s
 tts1n2

1
wi
s

2n

ˆi
s
n2
 n3

, (94)
i
s
 tts1n2

1
wi
s

2n

ˆi
s
n2
 n3

, (95)
where ts is given by Eq. 55 and
ˆi
s
wsi
s
/4 wsi
s

2/8. (96)
Insert Eqs. 94, 95, 49, and 50 into Eq. 56 and let n 3
to obtain
0 tH1pH0 wH0pH0
 ta1pH0efA
d
I
 ta1QI
ba1ecI . (97)
This equation can be simplified by noting that the Nernst–Planck equation,
Eq. 31, can be written in terms of dimensionless variables as
Js ts1pss ws(spss]. (98)
Substitution of this expression into Eq. 97 gives Eq. 57.
APPENDIX B: EQUILIBRIUM PROBABILITY FOR
PROTON OCCUPATION
This appendix shows how Eq. 70 is obtained from Eqs. 68 and 69. To
begin, consider the Boltzmann distribution of the density PH at equilib-
rium, Eq. 34 with s  H. From this equation, obtain an expression for
PH(H) in terms of PH(A
H) by eliminating KH. Integrate the result to find

A
H
A
H
PHH dH
 PHA
H
expWHA
H


A
H
A
H
expWHH dH.
(99)
Next, consider the second integral in Eq. 69. Use the Boltzmann distribu-
tion at equilibrium to solve for Pd(d) in terms of Pd(C
d ) and integrate to
obtain

C
d
C
d
Pdd dd
 PdC
d
expWdC
d


C
d
C
d
expWdd dd. (100)
We now consider QI
b. Start with Eq. 59 and substitute the middle
expression for a from Eq. 65. Transform the exponent using, from Eqs. 6
and 14,
f B
d
I 
d
C
d

d
B
d
, (101)
to obtain
QI
b
 PdC
d
expWdC
d
expdB
d


C
d
A
d
expdd dd. (102)
Finally, consider the probability QII
b appearing in Eq. 69. Obtain a new
expression by considering Eq. 61. The term Pd(C
d ) can be written in
terms of Pd(C
d ) by using the Boltzmann distribution for Pd, Eq. 34, and
eliminating Kd. Substitute the last expression for a from Eq. 65. The
exponent can be simplified using
f B
d
I 
d
C
d

d
B
d
 (103)
yielding
QII
b
 PdC
d
expWdC
d
expdB
d



A
d
C
d
expdd dd. (104)
26 Schumaker et al.
Biophysical Journal 80(1) 12–30

We now divide the right-hand side of Eq. 68 by the right-hand side of
Eq. 69 and substitute the results of Eqs. 99, 100, 102, and 104 to obtain
QFW
H


A
H
A
H expWHd

A
H
A
H expWHd
, (105)
where  is given by Eq. 71 and
 PdC
d
PHA
H

1expWdC
d
WHA
H
.
(106)
 can be transformed by substituting from Eq. 62 evaluated at equilibrium.
Further, use the definitions associated with the potentials s and s to find
WdC
d
WHA
H
C
d
A
H

I
2

f H  f C
d
I
2
,
(107)
leading to
 
H/dC•/CIexpC
d
A
H
I  . (108)
Simplifying the exponent using Eq. 44 finally yields Eq. 70.
APPENDIX C: GENERAL SOLUTION OF THE
FRAMEWORK MODEL
In this section, we find the general solution for the framework model
proton current. We will obtain a system of eight equations for eight
unknowns. Calculations are somewhat less cumbersome when made in
terms of the dimensionless variables defined in Appendix A. Further, many
equations can be written in parallel for the proton and defect species. We
will use the superscript s to denote either protons, s  H, or defects, s 
d. In this notation, the eight unknowns are the densities at either end of the
proton and defect diffusion intervals, namely ps(0) and ps(1), the proba-
bilities QI
b and QII
b of the defect boundary regions, and the proton and defect
fluxes, Js.
The first equation is simply JH  Jd, which must hold at steady state.
This is immediately eliminated by introducing a common variable for the
flux, J  Js. The boundary conditions, Eqs. 59 and 61–63, give us four
more equations. Written in terms of dimensionless variables, these are
Q I
b
 pd1a expf B
d
I, (109)
Q II
b
 pd0a expf B
d
I, (110)
pH0expf A
d
IJta  pd1expf B
d
I  cI ,
(111)
pH1expf A
d
I Jta  pd0expf B
d
I  cII , (112)
The sixth and seventh equations are obtained by integrating the Nernst–
Planck equation, Eq. 31, which, in terms of dimensionless variables, is Eq.
98. Multiply by an integrating factor and integrate to get
Jsts

0

s
expws  d  pssexpwss
 ps0expws0. (113)
Note that ws()  s() 
s(), where the
s() are given by Eqs. 92
and 93. It is convenient to define the integrals
hs

0

exps  f s2  1
I/2 d , (114)
where we use fd  fC
d and denote hs  hs(1). Then evaluate Eq. 113 at s 
1 and divide by a common factor to get
Jstshs  ps1expf s
I/2 ps0expf s
I/2. (115)
This is the form of the integrated Nernst–Planck equation to be used for the
general solution.
The eighth equation is obtained from the normalization condition,
which we write in terms of dimensionless variables

0
1
pHd 

0
1
pdd  QI
b
 QII
b
 1. (116)
We proceed to express the integrals in terms of J and ps(0). Solve Eq. 113
for ps(s) and integrate to get

0
1
psd Jsts

0
1
expws

0

expws d d
 ps0expws0

0
1
expwsd.
(117)
The following difference of energies will be useful,
ws  ws s 
f s2  1
I
2
 
s


fs2  1
I
2
. (118)
Then, the second term on the right-hand side of Eq. 117 can be rewritten
using
expws0

0
1
expws d  expf s
I/2g s, (119)
where we define the integral
gs 

0
1
exps f s2  1
I/2d. (120)
Further, define
is 

0
1
expws

0

expws d d


0
1
exps f s2  1
I/2hsd, (121)
Framework Model for H Conduction 27
Biophysical Journal 80(1) 12–30

where Eq. 118 and the definition, Eq. 114, are used to obtain the second
expression. Use this result and Eq. 119 to rewrite Eq. 117 as

0
1
psd Jstsis  ps0expf s
I/2gs. (122)
Substituting this expression and the boundary conditions, Eqs. 109 and
110, into the normalization condition, Eq. 116, we have the final form for
the eighth equation,
JtHiH  pH0expfH
I/2gHJtdid  pd0expf C
d
I/2gd
 pd1aexpf B
d
I pd0aexpf B
d
I 1. (123)
We now have a system of five equations and five unknowns. The
equations include the boundary conditions, Eqs. 111 and 112, the inte-
grated Nernst–Planck Eqs. 115, and the normalization condition, Eq. 123.
The unknowns are the proton and defect densities at the ends of their
intervals, ps(0) and ps(1), and the common current J.
These equations may be solved by straightforward algebra. For exam-
ple, eliminate pH(0) and pH(1) by substituting from Eqs. 111 and 112 into
Eq. 115 (for s  H) and Eq. 123. Then solve for pd(0) in Eq. 115 (for s 
d) and use this to eliminate pd(0) in the two remaining equations. This
leaves a system of two equations for the unknowns pd(1) and J. These can
be finally solved for J to obtain Eq. 74. In these calculations, there are
several points where Eq. 16 may be used to simplify an exponent.
In the denominator of J, the coefficient FH is given by Eq. 75, where
F0
H
 hHgd  2a coshf B
d
 f C
d /2
Ie
I/2,
(124)
F1
H
 iH  hHgHe
I, (125)
F2
H
 iH. (126)
The coefficient Fd is given by Eq. 76, where
F1
d
 id  hdgd  hda expf B
d
 f C
d /2
I]exp
I, (127)
F2
d
 id  hdaexpf B
d
 f C
d /2
I, (128)
F3
d
 gHhdexp
I/2 . (129)
Finally, the coefficient Fa is given by Eq. 77, where
F0
a
 2exp
I/2 gd  2a coshf B
d
 f C
d /2
I
 coshf A
d
 f H/2
I , (130)
F 1
a
 gHexp1 f A
d
 f H/2
I , (131)
F 2
a
 g Hexpf A
d
 f H/2
I. (132)
APPENDIX D: CURRENT SYMMETRY
In this appendix, we describe how the coefficient symmetries, Eqs. 79–81,
are obtained. These guarantee the current symmetry, Eq. 78, which must
hold by the symmetry of the gramicidin channel potential of mean force.
We work with the dimensionless variables described in Appendix A.
The coefficients Fi
s depend on the integrals hs, gs, and is, defined by Eqs.
114, 120, and 121, respectively. For this appendix, we define these inte-
grals as functions of  and
as follows:
hs,


0

exps  f s2  1
/2d , (133)
gs


0
1
exps  f s2  1
/2d , (134)
is


0
1
exps f s2  1
/2hs,
d.
(135)
As a special case, note hs(
)  hs(1,
). In the following, we make use of
the symmetries of the potentials of mean force, s()  s(1  ), and
change variables,   1   or   1  , in the integrands. We find
immediately
hs
 hs
, (136)
gs
 gs
, (137)
By the same method, it follows that
hs, 
 hs1,
 hs1 ,
, (138)
which is used to obtain
is
 gs
hs
 is
. (139)
Using these results, it is straightforward to show that the coefficients Fi
s,
given by Eqs. 124–132, have the symmetries of Eqs. 79–81.
APPENDIX E: SATURATING CURRENT AT
HIGH VOLTAGE
In this section, we obtain the asymptotic formula for the framework model
current as
I 3
, Eq. 82. Dimensionless variables are again used, see
Appendix A. We must first obtain asymptotic expressions for certain
integrals. The following calculations are intuitively reasonable and can be
justified rigorously by use of Watson’s lemma (Keener, 1988).
The integral hs() is defined by Eq. 114, where hs  hs(1). As
I 3
,
the integrand becomes sharply peaked near  0, and it is therefore
reasonable to write
hs
exps0

0

expf s2  1
I/2d

f s
I1expf s
I/2, (140)
where  means (f s
I)
1exp(f s
I/2) is asymptotic to h
s() as
I3
, in the
sense of asymptotic expansions. In particular, this result holds for hs 
hs(1). The integral gs is defined by Eq. 120. By an argument similar to the
one just given, we have
gs
f s
I1expf s
I/2, (141)
as
I3
. The integral i
s is defined by Eq. 121, depending on hs() whose
28 Schumaker et al.
Biophysical Journal 80(1) 12–30

asymptotic behavior is given by Eq. 140. As
I 3
, we have
is
f s
I1expf s
I/2


0
1
exps f s2  1
I/2d

f s
I2expf s
I. (142)
Finally, we must consider, as a special case, the combination gshs  is.
From the definitions, Eqs. 114, 120, and 121, we have
gshs  is 

0
1
d


1
d exps  s
 f s  
I. (143)
Making the substitution    , this becomes
gshs  is 

0
1
d expf s
IF, (144)
where
F

0
1
d exps   s. (145)
The integrand of Eq. 144 is sharply peaked near   0, leading to
gshs  is
f s
I1 (146)
as
I 3
.
We are now ready to obtain the asymptotic form of the framework
model current J in the limit
I 3
. The exact expression for J is Eq. 74,
where the terms in the denominator have the form given in Eqs. 75–77. The
numerator of J grows like e
I as
I 3
. The coefficients Fi
s in the
denominator are given by Eqs. 124–132. Using the asymptotic expressions
for hs, gs, is, and gshs  is, we find that the only coefficient that grows as
fast as e
I is F0
a, whose expression is given by Eq. 130. We find
F0
a

ae
I. (147)
Then, from the expression for J, Eqs. 74–77, we obtain the result expressed
by Eq. 82.
APPENDIX F: EQUILIBRIUM CONDUCTANCE
In this section, we obtain the expression, Eq. 84, for the framework model
conductance at a symmetrical equilibrium. Begin by substituting the ex-
pression for the current J, Eq. 74, into formula 83 for conductance. Taking
the derivative and then evaluating the resultant expression at equilibrium,
find that
dJ
dI

Eq

cIeI
tHFHtdFdtaFaEq
. (148)
When FH and Fd are written out, using Eqs. 75–77 and Eqs. 124–132,
they contain terms proportional to cIe
I
 cII, which vanish at equilibrium.
Introducing the definitions
A
I tHhH  2ta coshf A
d
 f H/2
I , (149)
B
I gd  2a coshf B
d
 f C
d /2
I . (150)
we can then write
tHFH  tdFd  taF aEq  exp 
I/2A
IB
I
 gHA
I t dhdB
IcII  tdhdgHexp 
I/2cIcII .
(151)
In the case of a symmetrical equilibrium, cI  cII  c and
I  0, the
denominator can be written in the factored form
tHFH  tdFd  taFaEq  A0e
I  tdh0
dcB0 g0
Hec.
(152)
Note that B(0)  gd  2a, where gd is defined by Eq. 120 and a by Eq. 65.
When gd is written in dimensional form, gd  2a  g1
d, as defined by Eq.
88. Inserting the result of Eq. 152 into Eq. 148 then leads to Eq. 84.
Mark F. Schumaker thanks Ann Schumaker for carefully editing the
manuscript. He is supported by grant MCB 9630475 from the National
Science Foundation.
Re´gis Pome`s is supported in part by the US Department of Energy through
the Los Alamos National Laboratory-Directed Research and Development
grant for bioremediation.
Benoıˆt Roux 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.
Chandrasekhar, S. 1943. Stochastic problems in physics and astronomy.
Rev. Mod. Phys. 15:1–89.
Chen, D., J. Lear, and R. S. Eisenberg. 1997. Permeation through an open
channel: Poisson–Nernst–Planck theory of a synthetic ion channel. Bio-
phys. J. 72:97–116.
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. Molec. Biol. 46:
51–96.
Dani, J. A., and D. G. Levitt. 1990. Diffusion and kinetic approaches to
describe permeation in ionic channels. J. Theor. Biol. 146:289–301.
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.
Dieckmann, G. R., J. D. Lear, Q. Zhong, M. L. Klein, W. F. DeGrado, and
K. A. Sharp. 1999. Exploration of the structural features defining the
conduction properties of a synthetic ion channel. Biophys. J. 76:
618–630.
Eisenman, G., B. Enos, J. Hagglund, and J. Sandblom. 1980. Gramicidin as
an example of a single-filing ionic channel. Ann. NY Acad. Sci. 339:
8–20.
Hille, B. 1992. Ionic Channels of Excitable Membranes, 2nd ed. Sinauer
Assoc. Inc., Sunderland, MA.
Hille, B., and W. Schwarz. 1978. Potassium channels as multi-ion single-
file pores. J. Gen. Physiol. 72:409–442.
Jackson, J. D. 1975. Classical Electrodynamics, 2nd ed. John Wiley &
Sons, New York.
Framework Model for H Conduction 29
Biophysical Journal 80(1) 12–30

Jakobsson, E., and S. W. Chiu. 1987. Stochastic theory of ion movement
in channels with single-ion occupancy: application to sodium perme-
ation of gramicidin channels. Biophys. J. 52:33–45.
Jordan, P. C., R. J. Bacquet, J. A. McCammon, and P. Tran. 1989. How
electrolyte shielding influences the membrane potential in transmem-
brane ion channels. Biophys. J. 55:1041–1052.
Karlin, S., and H. M. Taylor. 1981. A Second Course in Stochastic Pro-
cesses. Academic Press, New York. 157–191.
Keener, J. P. 1988. Principles of Applied Mathematics. Addison-Wesley,
Redwood City, CA. 434–437.
Kurnikova, M. G., R. D. Coalson, P. Graf, and A. Nitzan. 1999. A lattice
relaxation algorithm for three dimensional Poisson–Nernst–Planck the-
ory with application to ion transport through the gramicidin A channel.
Biophys. J. 76:642–656.
La¨uger, P. 1973. Ion transport through pores: a rate theory analysis.
Biochim. Biophys. Acta. 311:423–441.
Levitt, D. G. 1986. Interpretation of biological ion channel flux data:
reaction rate versus continuum theory. Annu. Rev. Biophys. Biophys.
Chem. 15:29–57.
McGill, P., and M. F. Schumaker. 1996. Boundary conditions for single-
ion diffusion. Biophys. J. 71:1723–1742.
Phillips, L. R., C. D. Cole, R. J. Hendershot, M. Cotten, T. A. Cross, and
D. D. Busath. 1999. Noncontact dipole effects on channel permeation.
III. Anomalous proton conductance effects in gramicidin. Biophys. J.
77:2492–2501.
Pome`s, R., and B. Roux. 1996. Structure and dynamics of a proton wire: a
theoretical study of H translocation along the single-file chain in the
gramicidin A channel. Biophys. J. 71:19–39.
Pome`s, R., and B. Roux. 1997. Free energy profiles governing H con-
ductance in proton wires. Biophys. J. 72:A246.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. 1992.
Numerical Recipes in Fortran. 2nd ed. Cambridge University Press, New
York. 272–273.
Robinson, R. A., and H. Stokes. 1965. Electrolyte Solutions. Butterworths,
London. 45–48.
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.
Roux, B. 1997. Influence of the membrane potential on the free energy of
an intrinsic protein. Biophys. J. 73:2980–2989.
Roux, B. 1999. Statistical mechanical equilibrium theory of selective ion
channels. Biophys. J. 77:139–153.
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, and B. Roux. 2000. A combined molecular
dynamics and diffusion model of single proton conductance through
gramicidin. Biophys. J. 79:2840–2857.
Zauderer, E. 1989. Partial Differential Equations of Applied Mathematics.
2nd ed. Wiley-Interscience, New York. 1–45.
30 Schumaker et al.
Biophysical Journal 80(1) 12–30