GLOSSARY

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS SUMMARY APPENDIX A APPENDIX B 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 R denotes either side I or side II of the channel.

Latin symbols

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

View larger version (19K): [in this window] [in a new window]   FIGURE 1   (A) Orientation of water molecules. Water molecules tend to be oriented so that each shares one hydrogen bond with the channel and two with neighboring waters. The axial component of the water dipole moments may reverse about a water that shares two bonds with the channel. This 2D depiction of hydrogen bonds is stylized; no true asymmetry about the pore axis is implied. (B) Applied potential. The single-proton model assumes that the electrical potential, V, is a linear function of the coordinate parallel to the pore axis, z. Side I refers to the solution on the left, and side II refers to the solution on the right. VI is the applied electrical potential on side I. VII = 0 by convention. (C) Hypothetical single-proton conduction mechanism. The top row of cartoons shows a pore with an excess proton. The proton may enter the channel on side I at the upper left. In this cartoon, the orientation of the pore waters favors entrance on side I. As the proton passes through the pore, the hydrogen bond structure of surrounding waters is altered so that water oxygens remain oriented toward the center of excess charge. The bottom row shows a pore in the absence of an excess proton. Water dipole moments flip as a defect in the hydrogen bonding structure diffuses through the column of waters. (D) State diagram of the single-proton model. The top segment corresponds to the proton-occupied pore and is parameterized by the coordinate µH. The bottom segment corresponds to a defect traversing the central barrier in Fig. 2 B and is parameterized by µd  [µCd, µCd]. The endpoints of the bottom segment are the boundary states bI and bII. They correspond to a defect in the boundary regions µCd < |µd| < µAd.

Greek symbols

	Friction coefficient (Appendix A)
	Framework model free parameter (see above Eq. 21)
µs	Dipole moment reaction coordinate for species s (see Eqs. 3 and 8)
±µAs	Maximum extents of dipole moment intervals (see Eqs. 6 and 10)
µw	Dipole moment of a water molecule (Eq. 1)
s	Scaling factor between the orientation and dipole moments for species s
s	Potential of mean force for species s (see Fig. 2)
s	Orientation moment for species s (see above Eq. 1)
mind	Orientation moment coordinate of d potential minimum
±As	Maximum extent of orientation moment interval (see Fig. 2 and Eqs. 5 and 9)
±Bd	Effective electrical coordinates of boundary regions (see Fig. 2 B)
±Cd	Maximum extent of interior of defect interval (see Fig. 2 B)
w	Orientation moment of a water molecule (Eq. 1)
	Electrical potential energy of pore contents (Eq. 23)

View larger version (20K): [in this window] [in a new window]   FIGURE 2   Proton and defect potentials of mean force. All energies are given in units of kBT for T = 298 K. (A) The ordinate is the proton potential, H, in kBT, and the abscissa is the orientation moment of the pore contents, H. Dipole moments are obtained by scaling µH = H H. (B) The ordinate is the defect potential, d, in kBT, and the abscissa is the orientation moment, d. Protons may enter the channel when defects are in the boundary regions I and II. Shown in the figure are wide boundary regions defined by Ad = 8.1, Bd = 6.8, and Cd = 5.7 e0Å. The text also considers narrow boundary regions with Ad = 8.1, Bd = 7.44, and Cd = 6.9 e0Å. Boundary regions surround the interior of the defect interval, defined by Cd  d  Cd. Dipole moments are obtained by scaling µd = dd.

	INTRODUCTION

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS SUMMARY APPENDIX A APPENDIX B REFERENCES

Proton conduction through gramicidin has gained recent attention because 1) Single-channel currents are large and can be measured over a very wide range of concentrations. The observed conductance has a very prominent shoulder as a function of [H⁺] between pH 1 and pH 2 (Eisenman et al., 1980). This observation has been confirmed by Akeson and Deamer (1991), and a similar shoulder has recently been observed in the RR dioxolane gramicidin analog by Cukierman (2000). 2) Experimental evidence suggests that the proton permeation mechanism is qualitatively different from that for other cations (Busath et al., 1998; Phillips et al., 1999). Molecular dynamics simulations suggest that water reorientation may be involved in the rate-limiting steps (Pomès and Roux, 1997, and manuscript in preparation). 3) Proton conduction can be used to probe permeation by other cations (Heinemann and Sigworth, 1989). 4) Proton conduction through gramicidin provides an interesting comparison with other proton channels (DeCoursey and Cherney, 2000). Among the latter are channels that form parts of molecular machines that convert proton gradients into chemical energy by synthesizing ATP (for example, Boyer, 1997) or convert proton gradients into mechanical energy by turning flagellar rotors (for example, Elston and Oster, 1997).

This paper develops a model for proton conduction through gramicidin based on the molecular dynamics simulations of Pomès and Roux (1997 and manuscript in preparation). These simulations incorporate the detailed structure of the channel (Arseniev et al., 1985). Techniques developed to describe the dynamics of this relatively simple molecule may soon be applied to channels of broader biological importance. Nuclear magnetic resonance and x-ray crystallography have recently provided high-resolution structures for several ion channels (Cowan et al., 1992; Ketchem et al., 1993; Doyle et al., 1998; Chang et al., 1998).

A large gap separates the time scales of molecular motions described by molecular dynamics and the time scales required to form meaningful averages of ion permeation events comparable to experimentally observed conductances. In a companion paper (Schumaker et al., 2000b) we develop a framework model for ion permeation that bridges the gap. This probabilistic model is designed to combine the results of molecular dynamics with simple assumptions about the proton entrance and exit processes that were not addressed by the simulations. The framework model is based on two coupled Nernst-Planck equations, one describing diffusion of a single proton and the other describing diffusion of a defect in the hydrogen bond structure of pore waters. The interaction between proton permeation and the dynamics of the pore waters can be understood in terms of a state diagram analogous to those encountered in rate theory.

The present paper is chiefly concerned with incorporating the results from molecular dynamics into the framework model and making a detailed comparison with experiment. We discuss the proton and water reorientation potentials of mean force and deconvolute their velocity autocorrelation functions to obtain numerical values for diffusion coefficients. We also verify that a mathematical approximation made by the framework model, the lumped state approximation, is accurate when applied to the water reorientation potential of mean force. After incorporating the information from molecular dynamics, the model has three free parameters. These are concerned with proton entrance and exitprocesses not addressed by the simulations.

Values for these parameters are obtained from the conductance and conductance ratio data of Eisenman et al. (1980), who present an extensive set of measurements for pH  1.7, where we believe single proton conductance dominates. A reasonable value is chosen for one parameter, and the other two are optimized to fit the data. A good fit is obtained, and the resulting model is also compared with two I-V curves measured by Decker and Levitt (1988). We also discuss our model in light of the recent experimental results of Phillips et al. (1999) on proton conduction in gramicidin analogs.

	METHODS

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS SUMMARY APPENDIX A APPENDIX B REFERENCES

Dipole moment of the gramidicin pore contents

The pore of gramicidin is roughly cylindrical, and consequently an experimentally applied transmembrane potential gives rise to an approximately constant electric field within the pore (Jordan, 1982; Jordan et al., 1989; Roux, 1999). A sketch of the pore geometry is given in Fig. 1 A. When the electric field is constant, the component of the pore potential energy due to the applied field depends only on the net charge and the z component of the dipole moment of the pore contents (Schumaker et al., 2000b). In part for this reason, we use the z component of the dipole moment of the pore contents as the reaction coordinate for the single-proton model.

The reaction coordinate of the molecular dynamics simulation of Pomès and Roux 1997; manuscript in preparation) depends on the orientation moment of the polarizable and dissociable PM6 waters (Stillinger and David, 1978; Weber and Stillinger, 1982) occupying the pore. Fig. 2, A and B, shows the potential of mean force of the pore contents with and without an excess proton, respectively, as a function of the z component of the orientation moment. This section describes what the orientation moment is and how it is rescaled to obtain an estimated dipole moment of the water column. We refer to a channel pore without an excess proton as an empty pore.

The orientation moment of a single water molecule is the formal dipole moment obtained by assigning oxygen a formal charge of 2e₀ and hydrogen a formal charge of +1e₀. The OH separation of the model water is ~1 Å, and the H-O-H angle is ~110°. Thus the orientation moment of the PM6 water molecule is ~2 cos 55°  1.15e₀ Å. Based on the work of Duca and Jordan (1997; see Results), we use the constant value of 2.4 Debye = 0.5e₀ Å for the water dipole moment in both empty and occupied pores. We therefore introduce a scaling factor of ^d = 0.5/1.15 between the orientation moment and the estimated water dipole moment.

A water typically forms one hydrogen bond with a pore backbone carbonyl and additional hydrogen bonds with the neighboring waters on either side. This configuration is suggested in Fig. 1, A and C, although in a highly stylized representation. The water dipole moments are tilted from the z axis by an average angle of ~45°. We use the value ^w = (1.15e₀ Å) cos 45°  0.81 e₀Å for the z component of the water orientation moment. Then the z component of the water dipole moment, µ^w, is given by

(1)

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

(2)

(3)

Then the relationship between the dipole and orientation moments of an empty channel is simply

(4)

Let n_max be the number of pore waters. To define the single file based on simulations, one can either count only those water molecules that can make up to two hydrogen bonds with other water molecules, or count only waters that are confined near the central axis of the channel pore. n_max = 10 is a consistent choice according to both criteria (Pomès and Roux, 1996). We must have n_max  n^d  n_max. Let _A^d and µ_A^d be the maximum magnitudes of ^d and µ^d, respectively, corresponding to configurations with all waters aligned. We then have _A^d  d  _A^d and µ_A^d  µ^d  µ_A^d, where

(5)

(6)

and µ_A^d = ^d_A^d. With the value of ^w obtained above, we have _A^d  8.1 e₀Å. In Fig. 2 B this value corresponds to the inflection points exterior to the potential minima on either side of the central barrier. Although the range of ^d values formally extends to ±9.1 e₀Å, the most extreme values correspond to strained configurations that are probably of little physical importance.

When the pore is occupied by an excess proton, the dipole moment includes a contribution due to the excess charge calculated about the center of the pore, z = 0:

(7)

(8)

where n^H is the net number of pore waters with dipole moments oriented in the +z direction in the occupied channel. We assume n_max  n^H  n_max. If the pore occupies the interval L/2  z  L/2, it follows that _A^H  H  _A^H and µ_A^H  µ^H  µ_A^H, where

(9)

(10)

In the occupied channel, the idealized state of maximum dipole moment corresponds to the excess charge at the extreme coordinate z = L/2, with n_max pore water dipole moments aligned in the z direction (see the upper right cartoon of Fig. 1 C). In this way, the partial negative charges carried by water oxygens are closest to the center of positive charge.

Combining Eqs. 5 and 9 and Eqs. 6 and 10, we have

(11)

From the molecular dynamics, _A^H = 3.35 e₀Å; then using _A^d = 8.1 e₀Å, we obtain L = 22.9 Å, ~2 Å shorter than the physical length of the pore. This discrepancy may be explained, at least in part, by the procedure used to simulate the channel occupied by an excess proton. In these simulations, the pore is occupied by 10 PM6 water molecules, and there are also several TIP3P waters clustered about the channel entrances at either end. The TIP3P waters cannot form hydrogen bonds. The closest that the excess proton can approach the channel entrances corresponds to hydronium-like (H₃O⁺) configurations on the outermost PM6 waters. If the TIP3P waters were able to host hydrogen bonds, there would be O₂H₅⁺ ions between the PM6 and TIP3P waters. Therefore, the value of _A^H obtained from the simulations may be slightly low, and, as a result, the effective length of the simulated channel may be slightly shorter than the pore of gramicidin.

Finally, we will obtain a scaling factor relating the molecular dynamics reaction coordinate ^H and the pore dipole moment µ^H. We make the simplifying assumption that when the excess proton is located at coordinate z, all water dipole moments will be oriented away from it, with water oxygens aligned toward the excess positive charge. Then n^H will be a linear function of z. Put

(12)

so that when z = L/2 (the proton is at the entrance on side II), n^H = n_max. Substituting Eq. 12 into Eqs. 7 and 8, and then dividing, we have

(13)

where

(14)

Diffusion coefficients

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

(15)

where µ is the reaction coordinate (the z component of the dipole moment of the pore content), m is a reduced mass, K is a spring constant, M is the memory kernel (a generalized friction coefficient), and a dot denotes the derivative with respect to time. F(t) represents the influence of degrees of freedom orthogonal to the reaction coordinate. Assuming that these other degrees of freedom relax on much shorter time scales than the reaction coordinate, F(t) is modeled as a random process.

If Eq. 15 is multiplied by (0) and an ensemble average is taken, an equation is obtained for the velocity autocorrelation function of µ. The result is then Laplace transformed, and the Einstein relation is applied to find a formula for the Laplace transform of the diffusion coefficient (Crouzy et al., 1994; Woolf and Roux, 1994):

(16)

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

(17)

For simplicity, we have assumed that the proton and defect diffusion coefficients are constant, independent of the value of their reaction coordinates. In both cases, c(t) was computed from molecular dynamics simulations at 0.5-fs intervals for times in the range 0 < t < 0.5 ps. Fig. 3 A shows the normalized proton velocity autocorrelation function ^H(t)^H(0) for the proton reaction coordinate held near ^H = 2.99 e₀Å by an umbrella potential. Values of the proton orientation moment are scaled to the dipole moment by µ^H = ^HH, as described in the previous section. Fig. 4 A shows ^H as a function of s. Parabolic extrapolation of ^H from the interval 10 < s < 40 gives an estimate ^H = 3.52 × (^H)² (e₀ Å)² ps¹. The smooth behavior of ^H at high values of s becomes corrupted for s < 10 corresponding to times t > 0.1 ps. This may reflect other processes influencing the transport of protons on longer time scales, for example, interactions between proton motion and rearrangement of hydrogen bonds.

View larger version (28K): [in this window] [in a new window]   FIGURE 3   Velocity autocorrelation functions. (A) The normalized proton velocity autocorrelation function H(t)H(0) for the reaction coordinate held near H = 2.99 e0Å by an umbrella potential. The unnormalized amplitude is H (t)2 = 7346.6 (e0 Å)2 ps1. (B) The normalized defect velocity autocorrelation function for the reaction coordinate held near d = 0.104 e0Å by an umbrella potential. The unnormalized amplitude is d(t)2 = 5270.8 (e0 Å)2 ps1.

View larger version (12K): [in this window] [in a new window]   FIGURE 4   Estimates of diffusion coefficients. In each figure, the numerical estimate of the Laplace transform of the diffusion coefficient, , is shown as a function of the transform parameter s. (A) Values of H(s). Parabolic extrapolation to the ordinate gives the estimated proton diffusion coefficient of H = 3.52 × (H)2 (e0 Å)2 ps1. The Laplace transform is shown for H = 1. (B) Values of d(s). The estimated defect diffusion coefficient is d = 1.08 × (d)2 (e0 Å)2 ps1. The transform is shown for d = 1.

Similarly, the defect diffusion coefficient was estimated from the defect reaction coordinate held near ^d = 0.104 e₀Å by an umbrella potential. Values of the defect dipole moment are given by µ^d = ^dd. Fig. 3 B shows the normalized defect velocity autocorrelation function. Fig. 4 B shows ^d as a function of s. Parabolic extrapolation of ^d from the interval 4 < s < 10 gives an estimate ^d = 1.08 × (^d)² (e₀ Å)² ps¹. Values of ^d become corrupted for s < 3 corresponding to times t > 0.33 ps; this reflects the finite interval of times for which c(t) was computed.

The lumped state approximation

The framework model lumps the boundary regions of the defect potential of mean force, shown in Fig. 2 B, into discrete boundary states. These are represented by the filled endpoints of the bottom segment of the state diagram shown by Fig. 1 D. The weights of these boundary states are carefully adjusted to give the integral over the Boltzmann distribution under the conditions of a symmetrical equilibrium (Schumaker et al., 2000b). This lumped state approximation (LSA) makes it possible to solve the framework model analytically. However, the boundary states act electrically like points with electrical coordinates ±_B^d, where _B^d is a free parameter of the framework model. The value of _B^d is chosen to give the best possible representation of the mean first passage time (MFPT) over the defect potential barrier as a function of the applied potential V_I.

To determine the best possible representation of the MFPT, consider a diffusion process that does not depend on the lumped state approximation. This latter process is defined on the interval [_A^d, _A^d] of the potential profile shown in Fig. 2 B. The process starts at ^d = _min^d = 6.5 e₀Å, a minimum of ^d, with a reflecting barrier at ^d = _A^d. The mean first passage time to ^d = _min^d is then calculated as described in Appendix B. Because the MFPT found in this way is the standard for comparison, we call it the exact MFPT. An analytical formula can also be obtained for LSA MFPTs, as will be described by Mapes and Schumaker (manuscript submitted for publication). These mean times correspond to a process that starts at the lumped state at ^d = _C^d. The exact mean first passage time to the lumped state at _C^d is then given. Values of the electrical coordinate for the lumped states, _B^d, are chosen so that LSA MFPTs approximate the exact values well over the range of applied potentials, 200 mV  V_I  200 mV.

Fig. 5 A compares the LSA MFPTs over the central barrier shown in Fig. 2 B with the exact MFPTs. The logarithm of the MFPTs is shown as a function of the applied potential V_I. When V_I > 0, the potential energy of the well centered at ^d = _min^d is increased relative to the energy of the well at _min^d, reducing the MFPT. To test the sensitivity of our results to the value of _C^d, MFPTs for the LSA are calculated for both wide and narrow boundary regions. The wide boundary regions are shown in Fig. 2 B and have _C^d = 5.7, _B^d = 6.8, and _A^d = 8.1 e₀Å. The narrow boundary regions have half the width of the wide regions, with _C^d = 6.9, _B^d = 7.44 e₀Å, and the same value of _A^d. Also shown in Fig. 5 A are MFPTs estimated by using Kramers' escape theory (e.g., Risken, 1989). The integrals in the Kramers formula are calculated numerically so that a robust estimate of the MFPT is obtained, even though the crest of the central barrier in Fig. 2 B has a complex shape; this is described in Appendix B. As shown by Fig. 5 A, MFPTs obtained by all four methods are in very good agreement over the full range of applied potentials considered.

View larger version (13K): [in this window] [in a new window]   FIGURE 5   Logarithm of the mean first passage time over the defect potential central barrier plotted as a function of the applied electrical potential. Squares and circles give estimates of the mean first passage time that apply to the potential defined in the full interval between ±Ad. Squares give the exact mean first passage between the potential minima, with reflecting boundary conditions at ±Ad. Circles give the Kramers approximation to the mean first passage time. Triangles give mean first passage times, using boundary states with Cd = 5.7 e0Å. Diamonds give times using boundary states with Cd = 6.9 e0Å. (A) Defect potential given by molecular dynamics. (B) Defect potential rescaled so that the height of the central barrier is 4kBT lower than that given by molecular dynamics.

The sensitivity analysis presented in the Results depends on the LSA remaining accurate for a range of barrier heights. In one case, the amplitude of the potential profile shown in Fig. 2 B is reduced by 4k_BT, leaving a barrier that is little more than 2k_BT high in the absence of an applied potential. Fig. 5 B shows how the LSA MFPTs then compare with the other MFPT values. Both sets of LSA MFPTs remain in good agreement with the exact values. However, estimates of the MFPT obtained from the Kramers formula diverge from the other values as V_I increases. This is not surprising, because Kramers' theory is an asymptotic approximation that is accurate when potential barriers are higher than ~4k_BT (Cooper et al., 1985, 1988). However, when V_I = 200 mV, the height of the potential barrier as seen by a diffuser beginning on side I is less than 1k_BT.

Fig. 5 shows that the logarithm of the mean first passage time is approximately a linear function of V_I. An interesting comparison can be made between linear least-squares fits to the points in the figure and the expected expression in the limit of high potential barriers. According to Kramers' theory, the rate k of barrier crossing in the latter limit has the form

(18)

where W^d is the height of the total defect potential barrier and A is a constant. The total potential, W^d, is the sum of the potential of mean force, ^d, and the electrical potential energy, ^d (Schumaker et al., 2000a). The mean first passage time corresponding to rate k is 1/k. The height, W^d, depends on the applied potential according to

(19)

where W₀^d is the height in the absence of an applied electrical potential and f is the fraction of the electrical potential drop between the maximum and minima of ^d. To be precise, put f = ^d_min^d/(e₀L), where _min^d = 6.5 e₀Å is the coordinate of the defect potential minimum and e₀L = 22.9 e₀Å is the total potential drop associated with the passage of a proton across the membrane (see Eq. 11). The predicted slope of log₁₀ (mean first passage time) as a function of V_I is then m, where

(20)

and where e is the base of the natural logarithms, e₀ is the elementary electrical charge, and the factor of 1000 is due to measuring V_I in mV.

Using the value of ^d described above (in Dipole Moment of the Gramicidin Pore Contents), we obtain, from Eq. 20, m = 0.00211. In comparison, linear least-squares fits through each of the four sets of values in Fig. 5 A give slopes in the range 0.00214  m  0.00216, close to the limiting value. However, linear least-squares fits through the exact and LSA sets in Fig. 5 B give 0.00189  m  0.00191. For the W₀^d  2k_BT potential barrier used to generate the latter figure, the LSA gives a softened dependence of mean first passage time on the applied field, which is similar to the exact result. In contrast, the slope obtained from the Kramers' estimates in Fig. 5 B gives m = 0.00209.

Optimization of model parameters

The framework model of proton conduction (Schumaker et al., 2000b) has three free parameters: , t^a, and _C^d.  depends on the absolute free energy difference between the proton and defect potentials of mean force, and t^a controls the rate of proton exit from the channel; together  and t^a control the rate of proton entrance and exit. _C^d determines the width of the boundary regions of the defect potential profile. These are defined as regions within which there may be a defect in the hydrogen bond structure of waters when a proton enters or leaves the pore (see Fig. 2 B).

Initial estimates of  and t^a are based on subjective estimates of the Michaelis constant, K_M, and the maximum conductance, G_max, associated with entry of the first proton into the gramicidin channel. K_M and G_max are assumed to be associated with the lower leg of rising conductance in the data of Eisenman et al. (1980), which is measured at symmetrical pH and a transmembrane potential of 50 mV (see Fig. 7 A). The estimated values are K_M = 0.005 M and G_max = 24 pS.

We use these estimates to determine values of X and Y in the formula for G₀, the conductance at symmetrical concentrations and zero transmembrane potential (Schumaker et al., 2000b):

(21)

In this formula, I is the proton current and V_I is the transmembrane potential. The subscript 0 refers to the derivative evaluated at the equilibrium state with symmetrical concentrations C_I = C_II = C and V_I = 0. The coefficients X, Y, and Z depend on the three free parameters in the following way: X = X(t^a, ), Y = Y (t^a, _C^d), and Z = Z (, _C^d).

Assuming a fixed value of _C^d, the estimates of X and Y determine initial values of t^a and . This is used as a starting point for local least-squares minimization of the model against the combined data in Fig. 7, A and B, for pH  1.7. t^a and  are varied to minimize the sum of unweighted differences between the log₁₀G₅₀ data (Fig. 7 A) and the model plus twice the sum of the unweighted differences between the conductance ratio data (Fig. 7 B) and the model. We consider two different values of the third parameter controlling the width of the boundary regions: _C^d = 5.7 e₀Å, giving wide boundary regions, and _C^d = 6.9 e₀Å, giving narrow boundary regions. These are used to generate the fits shown in Fig. 7, A and B. 

View larger version (38K): [in this window] [in a new window]   FIGURE 6   Single-proton model trajectories. Trajectories for A and B are obtained by simulating a random walk approximating the single-proton model for Cd = 5.7 and Ad = 9.1e0 Å, CI = CII = 0.01 M, and VI = 0. (A) A proton trajectory. (B) A defect trajectory. Occupation of the boundary states is indicated at Cd = ±5.7 e0Å. (C) Defect trajectories diffusing in the full interval of reaction coordinate defined by the molecular dynamics simulations.

View larger version (13K): [in this window] [in a new window]   FIGURE 7   Comparisons between the framework model and experimental data. Circles are the data of Eisenman et al. (1980). Curves show framework model fits to data. Fits are made optimizing the values of ta and , parameters that control the rate of proton entry and exit. Solid curves show fits using molecular dynamics estimates of potentials of mean forceand diffusion coefficients, and the wide boundary regions shown in Fig. 2 B; optimized values are ta = 23.5 ns and  = 3.60. Long-dashed curves show fits using a defect potential of mean force scaled to an amplitude 2kBT higher than that given by molecular dynamics; the proton potential of mean force and all diffusion coefficients are those given by molecular dynamics. Optimized values are ta = 16.5 ns and  = 3.53. Short-dashed curves show fits using molecular dynamics estimates of potentials and diffusion coefficients and narrow boundary regions. Optimized values are ta = 21.6 ns and  = 2.23. (A) Fits to conductance measured at 50 mV applied potential. (B) Fits to conductance ratios. The experimental data are conductances measured at 100 mV divided by estimated conductances at 0 mV. (C) Comparison between framework model fits to data of Eisenman et al. and the experimental I-V curve of Decker and Levitt (1988) at pH 3.75. (D) Comparison with Decker and Levitt I-V curve at pH 2.75.

	RESULTS

TOP ABSTRACT GLOSSARY INTRODUCTION METHODS RESULTS SUMMARY APPENDIX A APPENDIX B REFERENCES

Molecular dynamics simulations

Pomès and Roux (manuscript in preparation) have simulated the dynamics of proton conduction through gramicidin. The orientation of the channel pore in the membrane is sketched in Fig. 1 A. The packing of the pore contents in several configurations is represented in simplified form by Fig. 1 C. In one set of simulations the pore includes an excess proton; these are represented by the top row of the figure. The water oxygens tend to be oriented toward the center of excess charge. The resulting water dipole moments partially offset the contribution of the excess proton to the total pore dipole moment.

In a second set of simulations, there is no excess proton. These are represented by the bottom row of Fig. 1 C. It is energetically favorable for interior pore waters to share one hydrogen bond with the channel and two others with the adjacent water molecules. In this way they tend to have their dipole moments aligned with similar z components. However, pore waters near the channel entrances tend to be less ordered. As a result, there is often a localized disorder in the packing structure of the pore waters near the entrance, centered on a water molecule that shares two hydrogen bonds with the channel. This disordered region may diffuse from one end of the pore to the other, reversing the dipole moments of the pore waters.

Potentials of mean force

Diffusion coefficients

1

1

1

1

(22)

1

5

1

1

Single-proton conduction model

Hypothetical conduction mechanism

Applied transmembrane potential

(23)

Framework model