GLOSSARY

TOP ABSTRACT GLOSSARY INTRODUCTION CONSTRUCTION OF THE MODEL EQUILIBRIUM PROBABILITY FOR... GENERAL SOLUTION OF FRAMEWORK... PROPERTIES OF THE GENERAL... NUMERICAL SOLUTION DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F REFERENCES

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, â	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.
fAd, fBd, fCd = fd	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.
Qis	Probability that state si is occupied; see above Eq. 21.
QRb	Probability that bR is occupied; see above Eq. 32.
QFWH	Framework model probability that a proton occupies the channel; Eq. 68.
QSMH	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.

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.
is	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.
is	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, µis	Dipole moment reaction coordinate for species s. See Fig. 1 and below Eq. 17.
±µAs	Maximum extent of coordinate interval. See Fig. 1.
±µBd	Effective electrical coordinates of boundary regions. See Fig. 2 A and above Eq. 16.
±µCd	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, is	Dimensionless reaction coordinate for species s. See Eqs. 90 and 91.
Bd	See below Eq. 44.
Cd	Eq. 2.
AH	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

TOP ABSTRACT GLOSSARY INTRODUCTION CONSTRUCTION OF THE MODEL EQUILIBRIUM PROBABILITY FOR... GENERAL SOLUTION OF FRAMEWORK... PROPERTIES OF THE GENERAL... NUMERICAL SOLUTION DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F REFERENCES

This paper describes the construction of a framework model of single proton conduction through the ion channel gramicidin. A framework model is a kinetic model, designed to incorporate potentials of mean force and diffusion coefficients 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 coordinates of molecular dynamics simulations parameterize a simplified configuration space for the system being modeled. 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 channel, 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 example, one may easily describe channels whose occupancy is limited to a single ion (e.g., Lä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 Goldman-Hodgkin-Katz theory (reviewed by Hille, 1992; recently 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 corresponds 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 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).

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 distributed in time, corresponding to ions entering an empty channel at steady state.

The framework model we construct below is designed to incorporate the molecular dynamics results of Pomè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 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-contained. The boundary conditions we obtain allow the resulting model to be solved analytically.

The model is then solved for the special case of thermodynamic equilibrium, showing how the occupation probability 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 thermodynamic equilibrium. Finally, we describe the method of numerically simulating trajectories underlying the framework model. Schumaker et al. (2000) incorporate into the framework model the results of the molecular dynamics simulations of Pomès and Roux (manuscript in preparation). They then make a detailed comparison with experiment.

	CONSTRUCTION OF THE MODEL

TOP ABSTRACT GLOSSARY INTRODUCTION CONSTRUCTION OF THE MODEL EQUILIBRIUM PROBABILITY FOR... GENERAL SOLUTION OF FRAMEWORK... PROPERTIES OF THE GENERAL... NUMERICAL SOLUTION DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F REFERENCES

Figure 1 A shows a state diagram for the dynamics of proton permeation through gramicidin based on the results of Pomès and Roux (manuscript in preparation). The state diagram consists of two segments. The top segment corresponds 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. Pomè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 e₀L/2  µ^H  e₀L/2.

View larger version (13K): [in this window] [in a new window]   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 coordinate, µ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 Pomès and Roux (1997). Potentials are defined in the interval µAH < µH < µAH. (C) Defect PMF, d (dots), as calculated by Pomè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.

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 < e₀L/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 (Pomès and Roux, 1997, and manuscript in preparation). The simulated system included the channel, pore waters, and a few waters clustered outside each channel entrance. The membrane 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 potential 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 diffusion 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 conduction 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 between 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; compare with Fig. 1 C. The interior, corresponding to most of the central barrier, occupies the interval [µ_C^d, µ_C^d]. Surrounding this are boundary region I to the left and boundary 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.

View larger version (16K): [in this window] [in a new window]   FIGURE 2   State diagrams and random walks. (A) State diagram of hypothetical proton conduction mechanism; compare with Fig. 1 A. ±µAH and ±µAd 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 ±µCd. Boundary regions are the subintervals µCd < |µd| < µAd. In the framework model they have effective electrical coordinates ±µBd. The symbols fH and fAd, fBd, and fCd 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  . 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.

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 symmetrical about the midpoints µ^s = 0, s  {H, d}. We will assume that these potentials are exactly symmetrical. Further, there is an unknown energy difference between the 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,

(1)

(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, V_I, 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 = V_I/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),

(3)

where the electrical potential is given by the linear drop

(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,

(5)

(6)

where e₀ is the elementary electronic charge, and µ^s is the dipole moment reaction coordinate. Electrical potential energies corresponding to V_I = 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:

(7)

(8)

(9)

(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 e₀V_I. Adding Eqs. 7-10 and dividing by E gives

(11)

This equation relates the physical length of the pore with the maximum values of the proton and defect reaction coordinates. When the values µ_A^H and µ_A^d obtained from the molecular dynamics simulations of Pomès and Roux (manuscript 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 distances, analogous to those encountered in rate theory (e.g., Hille, 1992).

(12)

(13)

(14)

(15)

where the notation f^d is used in the appendices. These electrical distances are proportional to widths of subintervals 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

(16)

The total energy, W, of the pore contents is the sum of the intrinsic and applied potentials,

(17)

Random walk limit to a diffusion process

In this section, we construct the framework model for proton 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 integrals, the Nernst-Planck equations), which describe diffusion 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 diffusion 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 b_I and b_II 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 symmetrical 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 b_I and b_II 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. Transport 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 H_i, i  {1, 2, ... , n}, denotes a proton at coordinate µ_i^H = µ_A^H(2i/n  1), and state d_i, i  {1, 2, ... , n}, denotes a defect at coordinate µ_i^d = µ_C^d(2i/n  1). State s_i will refer to either H_i or d_i. The two additional states b_I and b_II will be taken to the boundary states of Fig. 2 B by the random walk construction. There are altogether 2n + 2 states in the random walk. We will define transition probabilities appropriately and take the limit n   to obtain the framework model.

Nernst-Planck equation for channel interior

i1

(18)

(19)

1

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

Entrance and exit transition probabilities

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

(44)

(45)

(46)