INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

The gramicidin A (gA) monomer is a pentadecapeptide consisting of alternating L and D amino acids. The conducting form of the channel is an N-terminal to N-terminal dimer composed of two identical ^6.5 helices. The channel is 25 Å long, with a central pore diameter of 4 Å (Arsenyev et al., 1990; Ketchem et al., 1997). The amino acid side chains extend radially outward, away from the helical backbone. These include four pairs of tryptophans at positions 9, 11, 13, and 15 positioned near the lipid-water interface. The tryptophan indole ring system has a dipole moment of ~2.1 Debye (Cotten et al., 1999b), comparable to that of a water molecule (Duca and Jordan, 1998). The orientations and dynamics of the indoles in the membrane and their dipole moment have been studied by solid-state NMR (Hu et al., 1993, 1995; Hu and Cross, 1995; Cotten et al., 1999b) and molecular modeling (Woolf and Roux, 1997; Dorigo et al., 1999; Anderson et al., 2001). The indole dipole moments contribute substantially to the electrical potential in the pore region. There is reasonable agreement (Hu and Cross, 1995) between the approximate magnitude of the dipole potential from the atomistic force field computations (Woolf and Roux, 1997; Dorigo et al., 1999; Anderson et al., 2001) and the effects of Trp-to-Phe changes on the Na⁺ conductance measured experimentally (Becker et al., 1991). However, the shape of the axial potential profile from the Trp side chains depends on the force field used. The results of Anderson et al. (2001), based upon an ab initio force field, are particularly simple. They find that the indoles of each tryptophan pair stabilize cations in the pore by ~0.6 kcal/mol. This change is approximately constant throughout the channel, independent of the spatial coordinate parallel to the pore axis, and presumably is extinguished by bulk electrolyte shielding just outside the channel.

An extensive set of gramicidin analogs have been developed and investigated by Andersen, Busath, Cukierman, Cross, Heitz, Koeppe, Woolley, and others (e.g., Andersen et al., 1998; Busath et al., 1998; Cotten et al., 1999b; Quigley et al., 2000; Jaikaran and Woolley, 1995; for reviews of earlier work, see Woolley and Wallace, 1992; Busath, 1993). In gramicidin B (gB), the tryptophans at position 11 are replaced by phenylalanine. The indole residues of the tryptophans are located outside the beta helix ~7 or 8 Å from the pore axis (Cotten et al., 1999b). The phenylalanine side chain is not expected to have a significant dipole moment. Therefore, it appears likely that differences between the conductance properties of gB and gA are mainly due to the change in the electrical potential in the pore region due to the decreased dipole moment of phenylalanine. Replacement of Trp by Phe provides a mechanism for tuning the electrostatic environment of the pore. In gramicidin M (gM), all four pairs of Trp are replaced by Phe (Heitz et al., 1982). The electrical potential of the gramicidin pore can also be modified by fluorinating the indole ring (Cotten et al., 1999b). In contrast to the effect of replacing Trp by Phe, 5-fluorination increases the side chain dipole moment according to both experiment (Andersen et al., 1998; Busath et al., 1998; Thompson et al., 2001) and computation (Anderson et al., 2001). Lipids also make an important contribution to the electrical potential within the pore interior. De Godoy and Cukierman (2001) report a recent study of this influence on proton conduction through dioxolane-linked gramicidin analogs. Measurements of proton conduction through gA in diphytanoylphosphatidylcholine (Rokitskaya et al., 2002) provide an intriguing comparison with the results in glycerolmonooleate reported here.

The models of proton conduction through gA, gB, and gM described in this work are based on the gA model of Schumaker et al. (2000; 2001). The system simulated by molecular dynamics included a representation of all of the atoms in gA, the pore waters, and a few waters clustered around each channel entrance. The calculations neglected interactions between the channel system and the surrounding lipid and most of the bulk water. Potentials of mean force and diffusion coefficients were calculated for two different sets of simulations. For both sets the reaction coordinate was essentially a scaled axial component of the dipole moment of the pore contents, calculated with respect to the center of the pore. One set of simulations modeled the occupied pore; an excess proton was included among 10 pore waters. In this case we designate the axial component of the dipole moment µ^H. The other set of simulations modeled the empty pore; the pore contents included only the 10 water molecules, and we designate the axial component of the pore dipole moment µ^d. The superscript d refers to defects in the hydrogen-bonded structure of the water chain that are thought to mediate reversal of the water dipole moments. The potential of mean force (PMF) from the simulations of the occupied pore is designated ^H and shown by the dots in Fig. 1 A. It has the form of a shallow potential well. In the absence of an applied electrical potential, the excess charge density is concentrated near µ^H = 0, close to the center of the pore. The PMF for the simulations of the empty pore is designated ^d and shown as the dots in Fig. 1 B. It is dominated by a central barrier centered at µ^d = 0, which we will refer to as the defect or water reorientation barrier. The probability density for µ^d is concentrated near the potential minima, on either side of the defect barrier. The model representation of applied transmembrane potentials ^H and ^d is indicated by the straight lines in Fig. 1, A and B. The assumption that an applied field is constant within the pore is reasonably accurate for the cylindrical geometry of the gramicidin pore (Jordan, 1982; Roux, 1999). In addition to the potentials shown, diffusion coefficients for the reaction coordinates µ^H and µ^d are obtained from their velocity autocorrelation functions (Schumaker et al., 2000).

View larger version (15K): [in this window] [in a new window]   FIGURE 1   (A) The intrinsic proton PMF H is calculated by molecular dynamics (Pomès and Roux, 1997) and the applied PMF H corresponds to a transmembrane potential of 150 mV. The reaction coordinate µH  [µ, µ] is the axial component of the dipole moment of the contents of the pore occupied by an excess proton. (B) The intrinsic defect PMF is d and the applied PMF is d. The solid curve shows the intrinsic PMF uniformly scaled so that its peak-to-peak amplitude is changed by  = 2.0 kcal/mol. The reaction coordinate µd  [µ, µ] is the axial component of the dipole moment of the empty pore (occupied only by waters). Boundary states correspond to the regions µ < |µd| < µ. Reorientation of pore waters is believed to be mediated by a packing defect that diffuses through the water column (Phillips et al., 1999; Pomès, 1999). (C) Model configuration space of single-proton conduction. The top segment represents the possible proton-occupied states and is parametrized by µH. The bottom segment represents the possible empty states and is parametrized by µd. Pairs of dashed lines delimit families of possible transitions between the proton and defect segments. They bound boundaryregions I and II on the defect segment. (D) State diagram of the single-proton model. Boundary regions I and II are lumped into boundary states bI and bII. The frequencies of transitions between the segments are indicated; R and R are rate constants, and CR is the proton concentration on side R  {I, II}.

These PMFs and diffusion coefficients were incorporated into a configuration space diffusion model of proton conduction through gA. The configuration space is analogous to the state diagram of a rate theory model. Our model does not make the mean field assumption, which is difficult to justify in the context of narrow pores (Hladky, 1999; Miller, 1999; Corry et al., 2000; Moy et al., 2000). The simplified configuration space is shown in Fig. 1 C. The top horizontal line segment represents states of the occupied pore and is parametrized by µ^H. The bottom horizontal line segment represents states of the empty pore and is parametrized by µ^d. The two pairs of dashed lines represent two possible sets of transitions between the occupied and empty segments of states; they define the boundary regions on either side of the empty segment. The cartoons at the four corners of the diagram represent extreme states at the ends of the line segments. The channel entrances are designated I and II, as indicated by the cartoon at the lower left. Consider an excess proton entering the channel on side I, at the upper left-hand corner of Fig. 1 C. The single-proton model assumes that entry can take place only when channel water dipoles are favorably aligned. As the center of excess charge diffuses through the pore, hydrogen bonds re-form so that the water dipole moments continue to point away from the excess charge. If the excess charge leaves the pore on side II, pore waters are left in an organized state. However, the waters are generally not perfectly aligned, and µ^d will be distributed in some region near the ^d potential minimum on side II. The µ^d reaction coordinate must then diffuse over the defect barrier so that the axial component of the pore water dipole moments reverse. Those waters will then be once again aligned to favor proton entry on side I.

This paper presents conductance data for gA, gB, and gM. The gramicidin analogs have dramatically different conduction properties. Consistent with the results of Eisenman et al. (1980), gA has two regimes of current increasing in proportion to concentration, separated by a shoulder between [H⁺] = 0.01 M and [H⁺] = 0.1 M. For gB, the shoulder is shifted by approximately half a pH unit in the direction of higher concentrations; gM has currents proportional to concentration from the lowest concentrations measured to [H⁺] > 1 M. No shoulder is visible in the gM current data. We explain these differences in the conduction properties of the gramicidin analogs as being due to differences in their electrostatic potentials in the pore region. Our model for the electrostatic potential difference between the analogs and gA is inspired by the electrostatic calculation of Anderson et al. (2001). We assume that the difference is constant, independent of the spatial coordinate parallel to the pore axis. We then find that the optimal value of this constant difference is very similar to the result of the electrostatic calculation. However, to obtain this result we have to decrease the amplitude of the water reorientation barrier in gA calculated by the molecular dynamics simulations of Pomès and Roux. A generalization of our result to Trp potentials that vary significantly in the pore interior is outlined in the Discussion.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

Experiment

Experimental setup and design was largely described in Phillips et al. (1999) with some exceptions. Single-channel conductance experiments were performed at various concentrations of HCl. Preliminary experiments were conducted using dilutions of 10 N HCl (Fisher Scientific, Fair Lawn, NJ). These experiments were later verified using exact concentrations of 0.9834 N HCl, 0.503 N HCl, 0.1933 N HCl, and 0.1006 N HCl obtained from Aldrich Chemical Company (Milwaukee, WI). Additional experiments were also collected at 0.01 N HCl, 0.02 N HCl, 0.05 N HCl, 0.002 N HCl, and 0.005 N HCl using 10- or 100-fold dilutions of the Aldrich HCl solutions. Experiments collected at 2 N HCl and 5 N HCl were diluted from the 10 N HCl Fisher solution. All dilutions were conducted using distilled water purified to >17.6 M-cm through a Barnstead NANOpure II system (VWR Scientific, San Francisco, CA). The combined data were further verified through direct comparison with the gA data of Eisenman et al. (1980).

Glyceryl monoolein (GMO, NuChek Prep, Elysian, MN) was dissolved in n-hexadecane (50 mg/ml) as described previously (Phillips et al., 1999). N-hexadecane was obtained from three separate sources, specifically from Spectrum Quality Products, Inc. (Gardens, CA, or New Brunswick, NJ), Fisher Scientific, and Aldrich Chemical Company, Inc. The lipid and solvent were used without further purification.

The gA was obtained from ICN Biomedicals (Aurora, OH); gM was prepared by solid-state synthesis using 9-fluroenylmethoxycarbonyl chemistry on an Applied Biosystems model 430a peptide synthesizer (Cotten et al., 1999a). The gM sample was recrystallized and/or repeatedly purified in 95% methanol and 5% water using an octadecylsilane column (Markham et al., 2001); gB was purified from gramicidin D (ICN Nutritional Biochemicals, Cleveland, OH) with an octadecylsilane column using 83% methanol/17% water (Koeppe and Weiss, 1981). Purified peptides were dissolved in HPLC grade methanol (Fisher Scientific) at ~10⁵ mg/ml. Approximately 10-20 pg of peptide were added to 2-3 ml salt baths to obtain the desired frequency of channel formation.

Lipid bilayer membranes were formed and channel conductance performance was measured as described previously (Phillips et al., 1999; Busath et al., 1998). Voltage was applied across the membrane using Ag-AgCl electrodes and a List Medical (Darmstadt, Germany) Patch Clamp L/M-EPC7, or a Warner BC-525C Bilayer Clamp (Warner Instruments, Hamden, CT). The data were collected using voltages between 25 mV and 150 mV at 25 mV increments. Experiments normally lasted between 15 and 40 min. Care was taken to avoid evaporative cooling artifacts. Currents were low-pass filtered at 30 Hz and sampled (100 per second) using a NI-DAQ Data Acquisition Board (National Instruments, Austin, TX) and IGOR Pro Version 3.12 software (Lake Oswego, OR) on a Macintosh computer. Channel conductance was analyzed using TAC X4.1 and TACfit V3.0 (Bruxton Corporation, Seattle, WA). Each experiment contained a minimum of 50 events, and at least three independent experiments were performed at each condition to verify results. Results were temperature-corrected to 23°C, as described by Phillips et al., 1999. Temperature deviations did not surpass 3°C. Data are available in tabular form at http://bioag.byu.edu/zoology/gramicidin. The model calculations described below were also performed assuming T = 23°C.

Model

Schumaker et al. (2000, 2001) developed a model of proton conduction through gA based on the molecular dynamics simulations of Pomès and Roux (1997). The latter authors simulated both the diffusion of an excess proton through the pore and the reorientation of water dipole moments in the absence of an excess proton. The reaction coordinates for both simulations can be directly related to the axial component of the dipole moments of the pore contents. The total change in dipole moment was equal to 22.9 e₀ Å. In principle, this result should be equal to an elementary charge passing through the length of the channel. Using the accepted value for the channel length (Urry, 1971; Wallace, 1990) gives ~25e_OÅ. A possible explanation for the discrepancy is that the proton entrance and exit process were not represented in the simulations.

When the original model with an effective length of 22.9 Å was adapted to a preliminary analysis of the gM data with the electrical distance associated with proton entrance optimized without constraint, a value was obtained that was greater than the total electrical width assigned to the proton entrance and exit steps (Schumaker, unpublished data). However, the analysis is made self-consistent when the channel length is increased to 25 Å and the additional electrical distance assigned to the proton entrance and exit steps. To do this, we begin by distinguishing between the length of the pore, L_P = 22.9 Å, and the total length of the channel, L_T = 25 Å. Appendix A shows how this distinction can be introduced into the formalism of the framework model (Schumaker et al., 2001) to give a proton conduction model with extended electrical distance.

With the boundary regions shown in Fig. 1 B, the original model has a reaction coordinate width of µ  µ = 0.55 e₀Å associated with proton entrance and exit on either side of the pore. Here, µ is an effective coordinate of the boundary region on side I (see Appendix A). Physically, this width corresponds to the increased order of pore waters in the presence of an excess proton at an entrance, as compared to the pore waters in the absence of an excess proton. The extended model has a width of f_E × e₀L_T = 1.60 e₀Å associated with proton entrance and exit at either side. The new electrical distance f_E includes both the effect of the increased orientation of pore waters and the contribution due to transport of a proton across the physical length of each entrance region.

The original single proton conduction model assigns the full electrical distance f to the exit rates _I and _II. This choice assumes that entry into gA is limited by diffusion to the channel mouth and not by a local barrier at the channel entrance. Diffusion into the channel mouth is not very sensitive to an applied electric field, except in the regime of low ionic strength and high applied fields where interfacial polarization becomes important (Andersen, 1983). Similarly, the present model for gA assigns the full electrical distance f_E to the exit rates. However, the models for gB and gM split f_E into two components:

(1)

where f_N is the component of the electrical distance assigned to the entrance step and f_X is the component assigned to the exit step. This allows for the possibility that proton entrance into the channel is limited by a local barrier near the entrance. Entrance rates should then be sensitive to an applied electric field.

As suggested by the calculations of Anderson et al. (2001) we assume that each pair of tryptophans in gA decreases the potential energy of a cation in the pore by a constant, independent of the ion's location along the pore axis. The Trp-11 pair is replaced by Phe in the gB dimer and all four pairs of Trps are replaced by Phe in the gM dimer. We model these replacements by a constant increase in the potential energy of an excess proton in the pore. Let ^HZ(µ^H) be the intrinsic component of the PMF of the proton-occupied gZ pore, where Z  {A, B, M} and µ^H is the proton reaction coordinate. The intrinsic component includes everything except the contribution of the applied transmembrane potential. The energy differences between the gramicidin analogs and gA are denoted by

(2)

where Z  {B, M} and ^B and ^M are independent of µ^H. Since we assume the electrical potential energy difference between the analogs and gA is constant, and the empty pore is electrically neutral as a whole, it also follows that

(3)

The increased potential energy of an excess proton in the pore of gB or gM has the effect of decreasing the probability of the occupied state as compared with that of gA under the same conditions. The changed probabilities of the occupied state must be reflected in changed values of the mean time before proton entrance into an empty pore and the mean time before proton exit from an occupied pore. Under the assumption that ^B and ^M are constant, these changed mean times must be due to changed rates of proton entrance and exit. To form a mathematical model, we decompose

(4)

where Z  {B, M}. ^NZ and ^XZ determine changes in the entrance and exit rates, respectively. This is described in more detail in Appendix B.

Optimization

Our measure of goodness of fit is _r, where the reduced chi-squared error is:

(5)

In this equation I is the kth experimental current value, I the corresponding single-proton model current value, and _k the experimental standard error of the mean. The number of degrees of freedom is  = N  M, where N is the number of data points and M is the number of model parameters optimized to fit the data (Press et al., 1992).

Fits of the single-proton model to the gA data are similar to those described by Schumaker et al. (2000). A fixed value of scales the defect PMF, as shown in Fig. 1 B. Data at four concentrations ([H⁺] = 2 mM, 5 mM, 10 mM, and 20 mM) and six voltages (V_I = 25 mV, 50 mV, 75 mV, 100 mV, 125 mV, and 150 mV) were used. Values of t^a and  were optimized, giving  = 24  2 = 22 degrees of freedom. Fits to the gB and gM data start with the results of the gA analysis. Values of the constant energy differences ^B and ^M, defined by Eqs. 2 and 3, are assumed; gB data at five concentrations (from [H⁺] = 2 mM through [H⁺] = 50 mM) and six voltages are used for the fit to the gB currents. Values of ^NB and f are optimized, giving  = 30  2 = 28 degrees of freedom. gM data at nine concentrations (from [H⁺] = 2 mM through [H⁺] = 1 M) and six voltages are used in the fit to the gM currents. Values of ^NM and f are optimized, giving  = 54  2 = 52 degrees of freedom. In this notation, f is the entrance electrical distance for analog Z  {B, M}.

Contour plots of _r as a function of ^NB and f and _r as a function of ^NM and f were examined in several cases, and multiple minima were never found in the physical domain of these parameters. In addition, the smooth variation of _r with parameters in Results/Sensitivity Analysis (see below) suggests that the optimized parameters are a smooth function of the assumed fixed values of ^B (or ^M) and . This suggests that jumps between multiple minima in the ^NZ × f plane (Z  {B, M}) did not occur in the construction of the sensitivity analysis.

The sensitivity analysis presented in Results determines the values of ^B or ^M and that were fixed in the optimization procedure described above. The optimization procedure was repeated on a grid of different values of ^B or ^M and . Intervals of  = 0.25 kcal/mol and ^B = 0.05 or ^M = 0.2 kcal/mol were chosen. The resulting grid of _r values was smoothly interpolated by a third-order polynomial in two variables to give a surface.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

Current measurements

Fig. 2 compares gA, gM, and gB single-channel current traces and conductance histograms under identical experimental conditions. The dominant peak in each histogram was designated the standard conductance for each peptide. With gramicidin M there seemed to be more heterogeneity in channel heights than with gA and gB. One consistent nonhomogeneity was a second peak ~30% of the standard peak current. A similar behavior was observed for this peptide previously in alkali metal cation solutions (Markham et al., 2001). For the purposes of this paper we assume that this second peak represents a minor conformational variant that can be ignored. According to our present data and the prior observation of Phillips et al. (1999), channel currents at 50 mV applied potential follow the sequence gA < gB < gM in 0.1 M HCl. This figure shows currents under similar conditions, except that the applied potential is 125 mV. The magnitudes of the currents follow gM ~ gA < gB.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   Comparison of representative gA, gM, and gB single-channel current traces and conductance histograms from individual experiments. All data were collected in GMO/hexadecane bilayers with 0.1 N HCl solution, 125 mV applied potential, at 21.5°C. At least 50 channels were observed for each experiment, and three or more experiments were performed at each condition. The dominant peak in each histogram was designated the standard conductance for each peptide. The top scale bars refer to the gA and gM current traces and the bottom scale bars to the gB trace.

Fig. 3 compares our gA data set with data taken from Eisenman et al. (1980). HCl and HNO₃ proton activities given by Eisenman et al. (1980) were scaled to concentration using the Mathematica package ActivityCoefficients.m developed by Arnd Roth (www.mathsource.com/Content/Applications/Chemistry/0206-198). Fig. 3 A shows gA channel current at 50 mV over a large range of log[H⁺]. The two sets of data are in good agreement. Both show a shoulder region between 0.01 and 0.1 M. Note that Eisenman et al. (1980) extended their measurements to very low proton concentrations using noise analysis, and found that currents in this regime are proportional to concentration (giving a slope of ~1 on a log-log scale). Fig. 3 B shows that the current-voltage relations shift from sublinear to superlinear near the onset of the shoulder just above 0.01 M in both data sets, peaking at 0.1 M and then declining toward 1.0 as saturation is reached. The maximum value of the conductance ratios is obtained at a somewhat lower concentration in our data than in those of Eisenman. In summary, the two data sets show the same qualitative features of current and conductance ratios.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Comparison of gA H+ currents at 50 mV (A) and chord conductance ratios (B) from previously published (stars, taken from Eisenman et al., 1980) and current (squares) experiments. Error bars in the current data represent ±1 SEM. In both experiments, bilayers were formed with GMO/hexadecane. The Eisenman et al. (1980) data were taken in the temperature range 22-26°C (see also Neher et al., 1978).

Fig. 4 A shows proton currents through gA measured over a wide range of symmetrical proton concentrations, from 2 mM to 5 M, and voltages from 25 mV to 150 mV. In the log-log representation there is a clear shoulder in the gA data with an inflection point at about log[H⁺] = 1.5. Above the shoulder there is a second regime with currents approximately proportional to concentration. Finally, for [H⁺] > 1 M, the rate of increase of conductance with concentration begins to level off.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Proton currents in symmetrical solutions as a function of [H+]. In the first three panels complete data sets for gA, gB, and gM are shown individually. Applied potentials of measured currents are indicated by symbols (bottom to top): 25 mV (diamonds), 50 mV (stars), 75 mV(squares), 100 mV (triangles), 125 mV (circles), and 150 mV (diamonds). Model currents at the respective voltages are indicated by solid curves. Fits by the single-proton model correspond to  = 2.0 kcal/mol; see Fig. 1 B. (A) Gramicidin A currents. Two parameters controlling the entrance and exit rates are optimized to achieve the fits. Optimized values are  = 1.959 kcal/mol and ta = 21.8 ns. (B) Gramicidin B currents. The parameters and B = 0.600 kcal/mol are fixed and the optimized values f  0 and NB = 0.046 kcal/mol are obtained. (C) Gramicidin M currents. The parameters and M = 3.090 kcal/mol are fixed and the optimized values f = 0.056 and NM = 0.634 kcal/mol are obtained. (D) Proton currents at the applied potential of 50 mV for gA (triangles), gB (stars), and gM (squares).

Fig. 4 B shows proton currents through gB measured over the same wide range of concentrations and voltages. The character of these data is generally similar to that of gA, but the shoulder is shifted by ~0.5 pH units, in the direction of higher concentrations. Below the shoulder there is clear convergence to a regime with currents proportional to concentration. The shoulder itself is especially well-defined at low voltages. There is little increase in current between [H⁺] = 100 mM and 200 mM at 25 mV. Above the shoulder there is again a second regime with currents approximately proportional to concentration. The transition from the shoulder to this second regime is very sharp, similar to observations made by Cukierman on RR dioxolane-linked gramicidin (2000). At concentrations above 1 M the rate of increases of conductance with concentration again levels off.

In comparison with the gA or gB currents, the gM currents shown in Fig. 4 C show reduced conductance at low concentrations and a greatly extended regime in which current is proportional to concentration. Similar to the other channels, the rate of increase of conductance with concentration levels off above 1 M. The apparent extension of the low concentration regime as tryptophan is replaced by phenylalanine admits a straightforward interpretation. If the low concentration regime (below the gA and gB shoulders) corresponds to conductance by a single-ion mechanism, and the high concentration regime (above the shoulders) corresponds to a multi-ion mechanism, the transition between these should depend on the electrical potential energy of a cation in the pore. According to all of the estimates (Woolf and Roux, 1997; Dorigo et al., 1999; Anderson et al., 2001) this potential energy should increase in the sequence gA  gB  gM. As the potential energy increases, the probability that the pore is occupied by a cation decreases, and the domain of single ion conductance should extend to higher concentrations.

Fig. 4 D compares the current-concentration relationships of gA, gB, and gM at 50 mV applied potential. At 0.002 M, gB (stars) starts out at a slightly lower conductance than gA (triangles), but then crosses over to a higher conductance at 0.005 M, reaching a shoulder at a higher concentration, and then eventually merging with the gA conductances. At 0.002 M the gM conductance is markedly lower than gA, and crosses over at ~20 mM. It never reaches a shoulder, but above 1 M the rise in conductance tapers off, as it does for the other two peptides. The decline in slope at high concentrations may reflect the same causal factors as with gA and gB. In particular, proton mobility in the bulk may be rate-limiting for [H⁺] > 2 M (Cukierman, 2000). However, the following discussion of conductance ratios suggests that the gM currents may incorporate a hidden shoulder.

To study the change in I-V curve shape with concentration we formed the conductance ratios g_V/g₂₅ for V = 50 mV, 100 mV, and 150 mV. In this expression, g_V = I_V/V is the chord conductance, where I_V is the current at applied voltage V. These ratios are shown in Fig. 5 for each analog and over the full range of concentrations measured. Conductance ratios less than one correspond to sublinear I-V curves and ratios >1 correspond to superlinear curves. Fig. 5 A shows the ratios for gA. Below the shoulder in Fig. 4 A the I-V curves are sublinear. There is a transition to superlinearity as [H⁺] increases from 5 mM to 50 mM, on the low concentration side of the shoulder. The conductance ratios are approximately linear between [H⁺] = 10 mM and 20 mM, near the inflection point in the log current versus log[H⁺] plot. The ratios attain their maximum values near [H⁺] = 50 mM and 100 mM (near the point of maximum curvature in the log current versus log[H⁺] plot) and then decline at higher concentrations.

View larger version (19K): [in this window] [in a new window]   FIGURE 5   Conductance ratios gV/g25 formed from proton conductances in Fig. 4. Symbols designate V: 50 mV (stars), 100 mV (triangles), and 150 mV (diamonds). Model conductance ratios are indicated by curves: 50 mV (solid), 100 mV (short dash), and 150 mV (long dash). (A) Gramidicin A. (B) Gramicidin B. (C) Gramicidin M. In C, thin lines indicate conductance ratios for a model fit assuming  = 1.25 kcal/mol.

Fig. 5 B shows the conductance ratios for gB. The general pattern is similar to gA, but is shifted toward higher concentrations by ~0.5 pH units. At low proton concentrations, I-V curves are sublinear; values of g₁₅₀/g₂₅ are as low as 0.6. The ratios reverse near [H⁺] = 50 mM. They attain a well-defined maximum at [H⁺] = 200 mM, at the same concentration as the very sharp transition to the regime above the shoulder in Fig. 4 B. Fig. 5 C shows the conductance ratios for gM. There is now a broad domain of sublinear I-V curves extending from the lowest concentrations measured to ~[H⁺] = 1 M. Compared with the corresponding regime for gB, ratios are not quite as low. g₁₅₀/g₂₅ remains >0.7. Conductance ratios are clearly increasing at [H⁺] = 500 mM and reverse near 1 M, consistent with the presence of a hidden shoulder.

Preliminary analysis of current measurements

In all three sets of data, I-V curves are sublinear at low concentrations and make a transition to superlinearity at a concentration that is fairly well-defined by the data. At sufficiently low concentrations, conductance should be limited by ion entry into the pore. If there is no significant local barrier to ion entry, then entry is limited by diffusion from the bulk to the channel mouth. This process is not expected to be strongly dependent on the transmembrane electrical potential, which would not extend far into the electrolyte solution. Thus, the sublinear I-V values at low concentrations can be readily understood. As concentration increases toward the shoulder region of the current data, the I-V curves make a transition to superlinearity. Evidently, a voltage-dependent step becomes rate-limiting. As we progress through the sequence gA  gB  gM, this transition takes place at higher concentrations, corresponding to higher currents and presumably higher proton entry rates. Thus, the data suggest that the voltage-dependent step becomes faster in the sequence gA  gB  gM. This line of reasoning can be continued to offer some support for the results of Anderson et al. (2001), who concluded that each tryptophan pair lowers the energy of a univalent cation in the pore of gramicidin by ~0.6 kcal/mol. If the onset to superlinearity corresponds to a transition from current limited by entry to current limited by a step taking place within the channel, then it should also correspond to rapidly increasing probability that the channel, as a whole, is occupied by a proton. This suggests that the proton concentration at the onset to superlinearity will be close to the proton dissociation constant. Let K_Z be the proton association constant and let E_Z be the difference between the energy of a proton in the channel and the energy of an empty channel in combination with a proton in bulk solution (Z = A, B or M). We should then have

(6)

and, as a consequence

(7)

where ^M = E_M  E_A, according to the notation introduced by Eq. 2. Approximating the proton dissociation constant with the proton concentration at the onset of superlinearity (these estimates are consistent with model calculations, see Eq. B14), we have from Fig. 5 that K  1 M and K  10 mM. From Eq. 7 we then have ^M 4.6k_BT  2.7 kcal/mol. The average cation energy contribution per tryptophan pair is then 0.7 kcal/mol, in rather good agreement with the results of Anderson et al. Fig. 5 B suggests that K_B  30 mM, with the result that ^B = E_B  E_A is about one-quarter of ^M. This would mean that the potential energy contribution of the Trp-11 pair is close to the average contribution per Trp pair, in further agreement with the results of Anderson et al.

We tentatively identify the domain of increasing conduction above the shoulders in Fig. 4, A and B with the onset of significant multi-proton conduction. This identification is consistent with our interpretation of the onset of superlinearity in Fig. 5, A and B as signifying rapidly increasing probability of single proton occupation of the gramicidin pore. The increase in occupation probability satisfies a pre-condition for multi-proton conduction, which requires ion entry into an already occupied channel. However, DeCoursey and Cherny (1999) have argued that the defect in the hydrogen-bonded chain that would necessarily exist between two excess protons in the pore would slow the rate of conduction below that observed experimentally. It would be very interesting to have molecular dynamics simulations of the multiply occupied pore, especially if the entrance and exit processes could be studied.

Single-proton model fits to current data

This section describes the fits of the extended single-proton model of gramicidin conduction to the current measurements of gA and its analogs shown in Fig. 4. The model is based on that described by Schumaker et al. (2000, 2001) but modified by increasing the effective length of the channel from 22.9 Å to 25 Å and assigning the difference to the channel entrance regions, as described in Methods and Appendix A. As suggested by the calculations of Anderson et al. (2001), we assume that each Trp pair makes a constant contribution to the electrical potential energy of a cation in the pore, independent of the spatial coordinate z parallel to the pore axis. The corresponding model variables are ^B and ^M (see Eq. 2); these are the positive energy differences between the proton-occupied states of the analogs and gA. The assumption that the contribution of the Trps is constant greatly simplifies the model, making the results independent of the detailed distribution of charge in the pore. The energy of interaction of the pore contents with the charge distribution of the indole ring system depends only on the net charge of the pore, which is +1 when the pore is occupied by an excess proton and 0 when the pore is only occupied by waters.

To obtain good fits to the gB and gM data it is necessary to reduce the height of the water reorientation barrier calculated by the molecular dynamics simulations of Pomès and Roux (1997). The original barrier is shown as the dots in Fig. 1 B. The present analysis assumes that the intrinsic potential ^d is uniformly scaled so that its maximum peak-to-peak amplitude is reduced by  =  2.0 kcal/mol. This scaled profile is shown as a solid curve in Fig. 1 B. The sensitivity analysis below gives an overview of how model fits depend on .

To achieve the fits to the gA data shown in Fig. 4 A, two parameters controlling the proton entrance and exit rates are optimized, holding the description of proton transport through the interior of the channel fixed. Proton entrance is assumed to be independent of applied potential, therefore all of the electrical distance associated with the channel entrances was assigned to the proton exit step. Only the data presented in Fig. 4 A between 2 mM and 20 mM were used in the fit, because currents too close to the shoulder might include a significant contribution due to multi-ion conduction. At each concentration, currents at all six voltages were used in the optimization. The agreement is good, except for some points at 125 mV and 150 mV. The variables that were optimized are  and t^a;  is equal to ^d(µ)  ^H(µ) plus a term that is proportional to the logarithm of the pore volume (Schumaker et al., 2001). t^a controls the proton exit rate (see Eqs. A15 and A16) and also determines the time scale of the mean time before exit (Appendix C). Together,  and t^a control the proton entrance rate (see Eqs. A13 and A14). Optimized values of  and t^a are given in the legend to Fig. 4.

Model currents in Fig. 4 A increase to a maximum value for [H⁺] between 100 mM and 500 mM, and then decline due to the clogging mechanism (Schumaker et al., 2000). To understand the clogging mechanism, consider an excess proton leaving the occupied pore on side II (Fig. 1 C). The waters of the empty pore must reverse their dipole moments before another proton can enter from side I. If the mean time required for proton entry on side II is shorter than the time required for water reorientation, a proton is most likely to re-enter the pore from side II. This frustrates the cycling around the diagram required for proton transport through the pore. Clogging reflects the finite rate of water reorientation and disappears as this rate is taken to infinity.

The fits to the gB and gM data shown in Fig. 4 use the gA analysis as a starting point. We assume a value for ^Z, Z  {B, M}, the constant difference between the energy of a cation in the pore of the analog and gA. Then two variables are optimized to make the fit to each analog. The first is the component of ^Z that alters the entry rate, ^NZ (see Eq. 4), and the second is the electrical distance associated with proton entrance, f.

The fit to the gB current data in Fig. 4 B is based on all of the data up through [H⁺] = 50 mM; ^B is the increase in energy of a cation in the pore of gB relative to gA. From the sensitivity analysis (presented below) this is fixed at the value ^B = 0.600 kcal/mol. This value is optimal for  = 2.0 kcal/mol under the additional constraint that f  0. It is in very close agreement with the calculations of Anderson et al. (2001). Optimized values of ^NB and f are given in the figure legend. The rate of entrance into gB is only ~8% less than the rate into gA. The electrical distance associated with entrance into gB is zero. The interpretation is that there is no significant activation barrier for proton entrance into gB and (consistently) the associated electrical distance is zero. One also obtains that the model gB exit rate is ~2.6 times the gA exit rate. In summary, our model results suggest that gB and gA share the property that there is no significant local barrier at the channel entrance to the entry of an excess proton into the empty channel. The difference between the conductance properties of the gB and gA models is almost entirely due to an increased rate of proton exit from the gB pore.

The fit to the gM current data in Fig. 4 C is based on all of the data up through [H⁺] = 1 M; ^M is the increase in the energy of a cation in the pore of gM, relative to gA, due to the replacement of all four Trp pairs with Phe. From the sensitivity analysis, the value ^M = 3.090 kcal/mol was fixed; this is optimal for  = 2.0 kcal/mol. Note that ^M  5^B. Optimized values of ^NM and f are given in the figure legend. The result for ^NM implies that the model gM entrance rate decreased by a factor of ~2.9, compared with the entrance rate into gA. Formally, this factor corresponds to an entrance barrier of ~1 kT. The associated reaction coordinate width assigned to the entrance step, f × e₀L_T = 1.41e₀Å, is only a little less than the total width, f_E × e₀L_T = 1.60 e₀Å, assigned to the entrance region. The value of ^XM obtained from Eq. 4 implies that the model gM exit rate is ~65 times the gA exit rate. To summarize, the model representation of gM includes substantially slowed entrance and much faster exit compared to gA or gB. Furthermore, the proton entrance process is voltage-dependent in gM, but not in gA or gB.

Comparison of model and experimental conductance ratios

Fig. 5 compares the observed conductance ratios with the proton conduction model. Values of the model parameters were obtained by optimizing the fit with the current data in Fig. 4 as described above. Model ratios were not directly optimized to fit observed ratios. The gA model ratios are compared with observed ratios in Fig. 5 A. The observed ratios are somewhat scattered; the process of forming g_V/g₂₅ leads to a significant standard error, especially at low concentrations. However, the model is in general agreement with the trend of the data below the proton concentration at which the conductance ratios reverse, ~20 mM. The low values of the model ratios in this regime (<1.0) are due to the rate-limiting entrance step, which does not depend on applied voltage. Model ratios continue to follow the observations to the maximum at [H⁺] = 100 mM, well above the highest concentration used in the optimization procedure. However, the model ratios fail to follow the subsequent decline in the observed ratios, as the model fails to follow the secondary rise in Fig. 4 A.

Fig. 5 B shows a particularly nice comparison between the observed and model gB current ratios below the concentration where the ratios reverse. Model ratios reverse for [H⁺] near 50 mM, in good agreement with the data. However, they do not quite match the very sharp increase in the observed ratios to a maximum at 200 mM. In Fig. 5 C, model gM ratios for  = 2.0 kcal/mol are the thick curves extending across the figure. These are in general agreement with the data for [H⁺]  500 mM. In particular, the model ratio g₁₅₀/g₂₅ saturates at a minimum level of ~0.8. This increase as compared to the model ratios g₁₅₀/g₂₅ for gA and gB is due to the voltage-dependent entrance rate. However, the experimental conductance ratios for gM reverse near [H⁺] = 1 M, while the model ratios reverse near [H⁺] = 2 M. Model ratios for  = 1.25 kcal/mol and the corresponding optimal value of ^M are shown as thin curves for [H⁺]  1 M. These reverse for a value [H⁺] > 5 M, clearly different from the experimental result.

Times required by the individual kinetic steps

Next, consider the times for the individual kinetic steps involved in the process of proton permeation from side I to side II according to the single-proton model. Suppose that a positive applied voltage V_I > 0 gives a net current of protons from side I to side II. Beginning at the upper left-hand corner of the state diagram in Fig. 1 D, the kinetic steps in the direction of the current are as follows. 1) The mean time t_X for a proton to exit the channel starting from just inside the channel entrance at side I. 2) The mean time t_R for a defect in an empty channel, starting at state b_II, to cross over the water reorientation barrier and reach state b_I. We calculate this with the possibility of the proton entering the channel on side II neglected. 3) The mean time t_N for a proton to enter the channel on side I, assuming that it remains in state b_I until that time. We denote the sum of these mean times by t_S. As the applied voltage on side I becomes large, the mean time for a proton to permeate the channel should approach t_S. In Fig. 6 we compare t_S with the current time, denoted t_C, which is simply the reciprocal of the current from side I to side II (the current in units of ions per nanosecond). An applied voltage of V_I = 150 mV was used to construct the figure. The times given correspond to the fits to the gA, gB, and gM data presented in Figs. 4 and 5. In particular, the water reorientation barrier shown in Fig. 1 B is uniformly scaled so that its amplitude is reduced by 2 kcal/mol.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   Times for individual kinetic steps in the model fit of Fig. 4. The applied potential is fixed at 150 mV. (A) Gramicidin A. From bottom to top and to the left of all crossover points, curves represent the water reorientation time tR (long dash), the exit time tX (solid), the entrance time tN (solid), the sum tS = tR + tX + tN (short dash), and the model current time tC (solid). Dots represent experimental current times. (B) Gramicidin B. (C) Gramicidin M. In the last panel tR  tX.

Fig. 6 A shows the different times for gA. At the lowest concentrations, near [H⁺] = 1 mM, the model current is limited by the entrance process. We may call this the entrance-limited regime. At moderate applied voltages t_C is substantially greater than t_S  t_N, because there is a significant probability that kinetic steps will take place in the direction opposed by the applied field. For example, an ion that has just entered the channel from side I may exit the channel on the same side. At very large applied voltages, t_C approaches t_N; and the current becomes more nearly diffusion-limited (results not shown). As concentration increases, t_N decreases in proportion to [H⁺]¹; t_N and t_Xcross near [H⁺] = 30 mM, at the low concentration side of the shoulder in Fig. 4 A. At slightly higher concentrations, the experimental values of the current time depart the model curve. This corresponds to the second regime of current proportional to concentration in Fig. 4 A. Throughout the concentration range, model reorientation times are much shorter than exit times.

The figure for gB is similar to that for gA, except that t_N and t_X cross near [H⁺] = 80 mM, at the low concentration side of the shoulder in Fig. 4 B. The gap between the exit and reorientation times is somewhat decreased as compared to that for gA, but these times are still well-separated. Fig. 6 C shows that the model representation for gM is entrance-limited for concentrations much below 1 M. The exit and reorientation time scales are now very nearly equal.

Spatial dependence of the mean time before exit

The mean time before exit t_X described in the previous section combines translocation across the pore interior with exit. However, the proton translocation process is itself very fast, due to the very high proton diffusion coefficient calculated by molecular dynamics (Schumaker et al., 2000) and the absence of a barrier in the proton PMF. The high proton translocation rate is apparent from Fig. 7, which shows the spatial dependence of the mean time before exit, (z), for gramicidin A and the analogs; z is the spatial coordinate parallel to the pore axis. Appendix C describes how (z) is calculated. Values of t_X correspond to evaluated at the left-hand endpoint of the curves in Fig. 7. The figure compares (z) for gA in the absence and presence of an applied potential, as well as comparing (z) for gA, gB, and gM in the presence of an applied potential. If the scale on these graphs were changed so that the base of the ordinate were  = 0, the graphs would all look nearly flat. The mean first passage time to escape the channel depends only slightly on the initial value of z within the channel. This reflects the different time scales to diffuse across the pore and to escape. According to the diffusion coefficient calculated from the molecular dynamics simulations (Schumaker et al., 2000), transport of the center of excess charge along a column of properly aligned water molecules is very fast. The mean time to leave the channel is much longer. As a consequence, model trajectories diffuse across the length of the pore many times before finally exiting. Fig. 6 A of Schumaker et al. (2000) gives a short segment of a typical trajectory.

View larger version (17K): [in this window] [in a new window]   FIGURE 7   Mean time before proton exit as a function of the coordinate z parallel t