Temperature Dependence of Gramicidin Channel Transport and

Article pubs.acs.org/JPCC

Temperature Dependence of Gramicidin Channel Transport and Structure Hyun Deok Song‡,§ and Thomas L. Beck*,†,‡ †

Department of Chemistry and ‡Department of Physics, University of Cincinnati, Cincinnati, Ohio 45221-0172, United States ABSTRACT: The gramicidin dipeptide is a small, cation-selective ion channel. Recent experiments have indicated that gramicidin can conduct ions at elevated temperatures. Since gramicidin is an efficient proton conductor, it is possible that this channel may have applications in fuel cell technology. In this study, we examine the temperature dependence of gramicidin A channel transport and structure with molecular dynamics simulations. In particular, the potentials of mean force (PMFs) for potassium ion motion through the channel are computed at five temperatures in the range 300−360 K. The channel displays a decrease in the free energy barrier height as the temperature increases. In addition, the enthalpic and entropic components of the free energy are computed, indicating a substantial enthalpy−entropy compensation and a positive entropy change when the ion enters the channel. The positive entropy change results from a reduction in fluctuations of ion interaction energies in the pore relative to those in the bulk solvent. The overall dimeric channel structure is maintained at 360 K for time scales up to 100 ns. In addition, higher temperatures affect the distributions of hydrogen bonds at the dimer interface, conformations of the N-terminal domain that may block the pore, channel bending angles, and distances between dimers. These findings may be related to the gating of the channel, although no complete dimer dissociation events were observed in the simulations.

■

INTRODUCTION Gramicidin is a hydrophobic dipeptide consisting of an alternating sequence of 15 D- and L-amino acids: (formyl-LX 1 -Gly 2 - L -Ala 3 - D -Leu 4 - L -Ala 5 - D -Val 6 - L -Val 7 - D -Val 8 - L-Trp 9 - D Leu10-L-Y11-D-Leu12- L-Trp13-D-Leu14-L-Trp15-ethanolamine). The dipeptide can act as an ion channel when embedded in a biological membrane. Gramicidin D is the natural mixture of gA (85%), gB (5%), and gC (10%).1 The three forms A−C are defined according to the residues X and Y as follows: gA1 (Val, Trp), gA2 (Ile, Trp), gB1 (Val, Phe), gC1 (Val, Tyr), gC2 (Ile, Tyr).2 There are multiple structural forms of the gramicidin A channel with two main types: a head-to-head single-stranded helical dimer3 (HD) known as the channel form and an intertwined double-stranded double helix4,5 (DH) known as the nonchannel form. Representative HD and DH structures are shown in Figure 1. Gramicidin is polymorphic in organic solvents, with the dominant form being the DH structure.6 DHs codissolved with lipid in organic solvent such as chloroform/methanol and ethanol are converted to HDs in the membrane environment after sonication and heat incubation, however, and this transition is temperaturedependent.7 High temperatures accelerate the transition,8 but the conversion is slowed at a high ratio of gramicidin to lipid.7 In simulation work, Siu and Böckmann9 found that the DH form displays an increase in ion conduction relative to the HD form due to a lower free energy barrier. It is believed, however, that the HD form is the active form in membranes,10−12 although some controversy has persisted about this issue.5,13−15 Recent NMR experiments have shown that the DH form can exist in membranes, but largely in aggregates;16 this indicates that, for channel-forming membranes relevant to gramicidin (such as DMPC), the predominant structure is the HD form. Thus, in this paper we examine the HD form in the simulations. © XXXX American Chemical Society

Extensive molecular dynamics simulations have been performed directed at the computation of potentials of mean force (PMFs) for ion motion through the HD form, with significant variation of the barrier height and PMF shape depending on the force field and free energy method used.9,17−25 Recent transport calculations for the HD form that include finite size and dielectric corrections result in conductance estimates in near-quantitative agreement with experiment.19 Another study found limited finite size effects on the PMF,26 and fully polarizable simulations indicate significant alterations in the ion PMF.25 Rather than apply computed PMF corrections as discussed above, in this study we scale the PMFs by a common factor that leads to agreement between computed conductances and experimental measurements at 300 K (and the PMFs at other temperatures are scaled by the same factor in the conductance calculations). Gramicidin selectively transports monovalent cations, with a conduction sequence H+ > NH4+ > Cs+ > K+ > Na+ > Li+.27 Recent experiments have indicated that the gramicidin channel can maintain its conductance properties at elevated temperatures near the boiling point of water.28−30 Cuppoletti et al. measured a proton conductance of 133 pS (in 1 M HCl) for the gramicidin channel in a 3:1 POPS:POPE bilayer at 298 K with a surrounding solution of 50 wt % glycerol; the conductance at 368 K was measured to be roughly twice the value at 298 K (unpublished results). These results may potentially have fuel cell implications, with the environmentally friendly gramicidin replacing Nafion as the membrane material. Received: June 6, 2012 Revised: February 6, 2013

A

dx.doi.org/10.1021/jp305557s | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

channel structure and transport. We end with a summary of our conclusions.

Fuel cells operate at temperatures similar to the highest temperature examined here.31 In the present paper, we model the transport of potassium ions through the gramicidin A channel, instead of protons, due to the difficulty of accurate representation of proton motion through the pore.32 We note that the room temperature conductance of H+ is around 30 times larger than for K+ through gramicidin in a DiPhPC membrane.33 The ratio of the ion diffusion constants in bulk water is 4.8.34 Clearly the mechanisms of proton and potassium ion transport through the channel are different, and thus the temperature dependence of the transport may also differ (see below). Here, we use the potassium ion as a simpler example for initial study, and we assume the same conductance scaling factor of 30 at the various temperatures in order to provide a rough comparison with the experiments discussed above. The scaling results in K+ conductance values smaller by roughly a factor of 6 than those in refs 33 and 35; the different membrane and solvent mixtures examined in the present study are likely the origin of these differences. Urry et al. 35 found that the gramicidin channel K + conductance increases as the temperature is increased from 288 to 320 K, with a break in the slope at roughly 300 K; the high temperature slope was found to be comparable to that for K+ diffusion in bulk water. By comparing with bulk transport data, they interpreted their results as rate-limited behavior at low temperatures (288−300 K) and diffusion-limited behavior at high temperatures (300−320 K); a reduction of the free energy barrier for ion permeation was discussed as a possible cause for the transition. In a more recent study, Chernyshev and Cukierman33 examined temperature-dependent conductivities for both H+ and K+ transport through the channel, and no break in the slope of the temperature-dependent conductance curve was observed for either of the ions (in the native gramicidin channel). In that study, the transport of both ions displays an activated form, but the enthalpy barrier for H+ transport (3.7 kcal/mol) is roughly one-half that for the K+ ion (7.2 kcal/mol). Also, negative entropies of activation for H+ transport were computed from the experimental data using Eyring rate theory. While these experimental studies provide interesting initial insights into the temperature dependence of the gramicidin channel, the underlying mechanisms of the temperature-dependent conductivity are not well understood. In the present paper, our goal is to examine the temperaturedependent behavior of the gramicidin A channel by computer simulation. As discussed above, we explore the structure and dynamics of the HD form, since a range of evidence suggests this is the active channel form in biological membranes. We have investigated the K+ ion PMFs at five evenly spaced temperatures over the range 300−360 K in a POPE bilayer to mimic the experimental system discussed above. (The CHARMM27 force field36 does not contain parameters for POPS lipids, so the mixed POPE/POPS membrane was not modeled here.) Most of our simulations also included 50 wt % glycerol as the external bath, also to mimic the experimental system. Some simulations were also performed for a pure water bath for comparison with previous modeling work.19 The paper is organized as follows. We next discuss the details of the simulations and the calculations of the ion PMFs, enthalpies and entropies, diffusion constants, and the maximum channel conductances. We then discuss the results of the simulations related to the temperature dependence of the

■

COMPUTATIONAL METHODS Simulation Protocol. The initial structures for the gramicidin A HD form was taken from the Protein Data Bank (1JNO37). The channel was embedded in 20 POPE lipid molecules that provide a single surrounding shell of bilayer molecules. The membrane/channel complex was then solvated with 698 TIP3P waters, 139 glycerols, and 4 KCl pairs with the VMD code,38 resulting in an ion concentration of roughly 0.3 M. The added glycerol is equivalent to a 50 wt % mixture in order to mimic the conditions in the experiments of Cuppoletti et al.28−30 Some simulations were also performed for a pure water bath using the same simulation protocol. The system was modeled with hexagonal periodic boundary conditions using the CHARMM27 force field,36 including the CHARMM27 lipid parameters (CMAP corrections were not applied39). A time step of 2 fs was used. After energy minimization by the conjugate gradient method, the system was heated to the desired temperature. This was followed by 5 ns equilibration and 4.5 ns runs for each umbrella sampling window (500 ps of equilibration following placement of the K+ ion and 4 ns of data collection).

Figure 1. Helical dimer (A) and double-helix (B) forms of gramicidin A (PDB IDs: 1JNO and 1AV2, respectively). NewCartoon expressions were generated with the VMD software.

The Langevin piston method40 was utilized for temperature control with damping coefficient 5 ps−1, and the pressure was maintained semi-isotropically (as was done in ref 9) at 1 atm with a modified Nosé−Hoover Langevin piston. The impact of the pressure coupling was examined by monitoring the average hexagon area of the simulation box over the course of a 10 ns simulation at 300 K. The area stabilized at a mean value of 775 Å2. From this value and the estimated area of the gramicidin channel (232 Å2),41 the calculated area per lipid is 54 Å2, close to the experimental value of 56 Å2.42 The computed area per lipid increased by about 6% when the temperature was increased to 360 K. A previous study43 concluded that applying either isotropic pressure or constant surface tension boundary conditions had little effect on the simulation results. The B



Article

particle-mesh Ewald method44 for long-ranged electrostatic calculations was employed. Hydrogen bond lengths were fixed with the SHAKE method.45 Nonbonded van der Waals interactions were smoothly decayed at 8.5 Å and truncated at 10.0 Å. The NAMD 2.6 simulation code46 was employed for all of the molecular dynamics simulations. During the simulations, the center of mass of the channel was held fixed with a strong harmonic constraint; this constraint only serves to center the channel and has no impact on the results. Also, as in ref 19, a weak harmonic constraint (5 kcal/ (mol Å2)) was applied to the z-component of the lipid center of mass. Position data were output every 2 fs for umbrella sampling calculations. Position and velocity data were saved every 10 fs for diffusion coefficient calculations. In order to check whether the system was fully equilibrated, the root mean square displacement (RMSD) from the experimental geometry was calculated. Once the RMSD stabilized, it was assumed that the system was in a basin of the free energy, which corresponds to quasi-equilibrium conditions.47 The average RMSDs for the HD form range from 0.70 Å at 300 K to 0.76 Å at 360 K. Potentials of Mean Force (PMFs). The cation permeation process through the gramicidin channel requires that the ion surmount a free energy barrier of several kBT. This barrier significantly decreases the probability of spontaneous ion transit, and this in turn makes the sampling difficult.48 This problem can be overcome by employing the umbrella sampling method, which introduces a biasing potential to improve sampling efficiency. The bias can subsequently be removed with the weighted histogram analysis method (WHAM).49−51 We chose the reaction coordinate as the z distance perpendicular to the membrane plane from the center of mass of the protein to the position of the potassium ion.19 Various methods can be used to define the reaction coordinate, which may slightly change the height of the free energy barrier. The absolute reaction rate is not affected by a change of the microscopic order parameter, however.19,52 If a water molecule was close to the center of the window inside the channel, the water was exchanged with a K+ ion in the bulk region. The exchanged K+ ion was dragged to the center of each window and then sampled with the biasing harmonic potential wi(z) = 1/2k(z − z0)2, where z0 is the center of the window and the spring constant, k, was chosen as 10 kcal/(mol Å2). This choice of the force constant yielded sufficient overlap of the probability distributions. In the bulk region, the K+ ion was simply pushed to the center of the window. A total of 81 windows from −20 to 20 Å, in increments of 0.5 Å, were generated. A flat-bottom harmonic potential with radius of 5 Å and a force constant of 10 kcal/ (mol Å2) was applied to the K+ ion in the bulk region (z values with magnitudes greater than 14 Å from the channel center of mass) to minimize lateral motion of the ion. In addition, by applying a spherical restraint with force constant 10 kcal/(mol Å2) and radius of 14 Å from the channel center of mass, other ions were not allowed to enter the channel during sampling, thus ensuring single ion occupancy at all times. The first 500 ps equilibration period of each window simulation was discarded, and the biased distribution of the ith window ⟨P(z)⟩(i) was calculated from the subsequent 4 ns simulation. The PMF, W(z) = Δμex(z), is the free energy (or excess chemical potential) profile along the reaction coordinate. The PMFs were obtained by unbiasing the umbrella sampling distributions using the WHAM method.49−51 A total of 800 bins with bin interval 0.05 Å were generated. The bootstrap

method in the WHAM code was utilized to estimate the errors in the computed PMFs. The WHAM calculations were iterated until the change of the free energy constants was less than 10−5 kcal/mol. Enthalpies and Entropies. The enthalpy and entropy profiles (Δhex(z) and Δsex(z), respectively) can be obtained from the temperature-dependent Δμex(z) (PMF) values as follows: ⎛ ∂β Δμex (z) ⎞ Δhex (z) = ⎜ ⎟ ∂β ⎝ ⎠P

(1)

where β = 1/kT and ⎛ ∂Δμex (z) ⎞ ⎟ Δs ex (z) = −⎜ ⎝ ∂T ⎠ P

(2)

The enthalpy and entropy profiles are computed directly from the temperature-dependent PMFs obtained from the simulations using a simple finite-difference approximation.53 In order to gain some insight into the origins of the enthalpy and entropy profiles discussed below, we further examine the statistical mechanical expression for the excess chemical potential:54 βμex = −ln⟨exp( −βε)⟩0

(3)

Here ε = UN+X − UN − UX is the interaction energy of the ion X with all of the surrounding molecules; UN+X is the total potential energy of the solvent plus ion system, UN is the potential energy of the solvent, and UX is the potential energy of the isolated ion, which is 0 for monatomic ions. The sampling in eq 3 is conducted with the ion absent from the system, indicated by the 0 subscript. Alternatively, the excess chemical potential can be written as βμex = ln⟨exp(βε)⟩

(4)

in which the sampling is conducted with the ion present. We use this coupled form for the analysis presented below. Taking the appropriate temperature derivative of the excess chemical potential, the exact solvation enthalpy is hex = ⟨ε⟩ + USR + Pv

(5)

where USR = ⟨UN⟩ − ⟨UN⟩0 is the change in the average potential energy of the surroundings upon adding the ion (solvent reorganization energy in bulk solution or in the protein sites), and the last Pv term involving the partial molar volume is negligible for this study. The solvation entropy is then s ex = − k ln⟨exp[β(ε − ⟨ε⟩]⟩ + USR /T

(6)

The partial molar entropy is thus the sum of a fluctuation term and the solvent reorganization contribution. Note that the USR term cancels exactly when constructing the free energy at a given temperature and thus does not contribute to the free energy profile; nevertheless, this term is essential for the separate enthalpy and entropy contributions and thus the temperature dependence of the free energy. While accurate methods exist for computing ion free energies,55 Gaussian models have provided helpful insights into ion selectivity in channels56 and solvation phenomena more generally.57 A cumulant expansion of eq 4 yields C


The Journal of Physical Chemistry C μex ≈ ⟨ε⟩ +

β 2 ⟨δε ⟩ 2

Article

Laplace transform of the memory function, we can obtain the following form for the diffusion constant as a function of s:60

(7)

where δε = ε − ⟨ε⟩. The corresponding Gaussian estimate for the entropy is s ex ≈ −k

−Ĉ(s)⟨vz(0)2 ⟩⟨Δz(0)2 ⟩

D̂ (s) =

(

Ĉ(s) ⟨Δz(0)2 ⟩s +

2

β ⟨δε 2⟩ + USR /T 2

(8)

2

D = lim D̂ (s)

(11)

(12)

s→0

+

Computed diffusion constants for the K ion in bulk water are 0.19 Å2/ps (MSD), 0.17 Å2/ps (VACF), and 0.20 Å2/ps (GLEHO). Our results are very close to previous results.60 The Maximum Conductance. The one-dimensional steady-state flux (J) of an ion in the Nernst−Planck continuum theory is given by61,62 J = −D(z)

(9)

dP(z) D(z) dWtot(z) − P(z) dz kBT dz

(13)

where Wtot is the total PMF that is the sum of the equilibrium PMF W(z) and the interactions due to coupling between atomic charges and the transmembrane potential ϕmp.19,63 This approach assumes that the external field does not appreciably alter the underlying PMF; that is, we are making a small-field approximation.62 It is also assumed that at most one ion occupies the channel. A constant field assumption was made, equivalent to assuming a linear applied voltage variation inside the channel; this assumption is equivalent to assuming a negligible charge rearrangement in the bath with changing voltage. The conductance G0 is

while the Gaussian estimate of the entropic part of the free energy change is β T Δs ex ≈ − Δ⟨δε 2⟩ + ΔUSR 2

) − ⟨Δz(0) ⟩

where C(t) is the normalized velocity autocorrelation function, and the Laplace transformation of C(t) is expressed as Ĉ (s) = −st ∫∞ 0 dt e C(t). The diffusion coefficient can be computed by taking the limit as s goes to zero (see below for further discussion of the limit)60

For ion hydration, the excess entropies are typically remarkably small in magnitude,58 comparable to those for rare gases. On the other hand, in a Gaussian model the hydration free energy is roughly −β⟨δε2⟩/2 ≈ ⟨ε⟩/2. This suggests that the solvent reorganization term USR is on the order of −⟨ε⟩/2, or close to −μex, and is thus positive with a substantial magnitude. This also gives a qualitative explanation of why hydration enthalpies are similar to the hydration free energies in water. The above crude estimate of the solvent reorganization energy is given only to indicate its large positive magnitude in a given environment, since it assumes a solvation entropy of 0, and we do not use this approach for estimating changes in the solvent reorganization energy below. The resulting enthalpy change for the ion moving from bulk into the pore is thus Δhex = Δ⟨ε⟩ + ΔUSR

⟨vz(0)2 ⟩ s

(10)

In the calculations presented below, finite-size electrostatic selfenergy corrections in periodic boundaries54 are not applied since they do not vary with ion position, and only changes between bulk and the channel are considered here. Dielectric corrections (that mimic the dielectric environment of the membrane), on the other hand, may lead to a significant reduction of the barrier height,19 but as mentioned above, rather than computing those corrections, we rather scale the PMF to obtain agreement with the experimental conductance at 300 K. These points will be discussed below. Calculations of the interaction energies and their fluctuations were performed at 300, 330, and 360 K during 10 ns simulations. The ion was constrained to a particular z location with a force constant of 10 kcal/(mol Å2). The energies were computed at 0.1 ps intervals, resulting in 105 interaction energies that were averaged. This data set allows for accurate estimates of mean interaction energies (to within less than 0.1 kcal/mol) and solid estimates of the fluctuation contributions (to within roughly 2 kcal/mol) but does not allow for accurate calculations of higher terms in the cumulant expansion in the pore. (Consistent results were observed for the third-order term in the bulk solution on both sides of the channel, however.) Diffusion Constants. In order to predict channel conductance, the spatially dependent diffusion constant is required. To start, we calculated the 1D diffusion coefficient of the K+ ion in 216 TIP3P waters by three different methods: mean-squared displacement (MSD), velocity autocorrelation function (VACF), and generalized Langevin equation59 for the harmonic oscillator (GLE-HO). The GLE-HO diffusion constant is computed as follows. By using the Einstein relation, the equipartition theorem, and the

G0 =

=

I ϕmp

(14)

zeJ RTϕmp ̅ /zF

(15)

where R is the gas constant, F is the Faraday constant, z is the valence, and ϕ̅ mp is the dimensionless applied potential. With the assumption of equal ion concentration on both sides of channel, eq 13 can be integrated, and J can be inserted in eq 15. Experimental data display linear I−V curves for the gramicidin channel.64 Under the small-field assumption, the maximum conductance is given by61,62 Gmax =

e2 ⟨D(z)−1e+W (z)/ kBT ⟩−1⟨e−W (z)/ kBT ⟩−1 kBTL2

(16)

where L is the effective length of the channel and the brackets imply a spatial average. The value L = 28 Å was chosen because deviation of the PMF from the bulk value started at roughly ±14 Å (see below). The leading term of eq 16 is the first bracket. In the limits of nonvarying diffusion constants in the pore and high PMF barriers, the conductance is linearly proportional to the pore diffusion coefficient, and it depends exponentially on the PMF barrier height. D



■

Article

RESULTS AND DISCUSSION PMFs. Figure 2 displays the computed PMFs for the gramicidin A HD form (water/glycerol and water baths). The

interested in the temperature dependence of the ion transport for a given model. The number of lipid bilayer molecules and the force field parameters can affect the height of the PMF barrier.9,18−24 In ref 19, when 96 DMPC lipids were inserted around the gramicidin channel, the PMF barrier decreased by about 1 kcal/ mol relative to when 20 lipid molecules were used. When the GROMOS87 force field was employed, the resulting PMF barrier increased by 4 kcal/mol relative to the barrier obtained using the CHARMM27 force field.19 The more recent GROMOS96 force field yields better agreement with the CHARMM27 parameter set.9 In order to compare with previous simulation results, we ran simulations at 330 K for the pure water bath, and our results are compared with other research results in Table 1. We find that the present results (for Table 1. Comparison of the HD PMF (Pure Water Bath) with Previous Results (kcal/mol)a

Figure 2. PMF profiles of the gramicidin A HD channel (A, B). PMFs (A) were symmetrized (B) as discussed in the text. Most simulations were run in a water/glycerol (w/g) bath at 300, 315, 330, 345, and 360 K. One extra simulation at 330 K was run in a water (w) bath in order to compare with previous studies. The error bars calculated using the Monte Carlo bootstrap method are not shown. The estimated sizes of error bars for asymmetric and symmetric PMFs are about 0.2 and 0.1 kcal/mol, respectively.

study

force field

no. of lipids

barrier

well depth

total barrier

this study (330 K, 2 ns) this study (330 K, 4 ns) ref 19 (330 K) ref 19 (310 K) ref 22 (298 K)

CHARMM27

20

6.5

−2.4

8.9

CHARMM27

20

6.0

−2.3

8.3

CHARMM27 GROMOS96 CHARMM27

20 124 128

9.6 9.3 9.1

−1.4 0 −2.9

11.0 9.3 12.0

a

Note that the barrier height decreases somewhat with increasing sampling time for each umbrella sampling window, as was observed in ref 19. The barrier heights reported here are for the PMF value at z = 0.

PMFs calculated directly from MD simulations are asymmetric with a deviation of roughly 1 kcal/mol on the two sides (Figure 2A). The main reasons for this asymmetry are thought to be incomplete sampling, fluctuation of the center of mass of the protein, membrane undulations (and local lipid fluctuations), and asymmetric structures of gramicidin around the center of mass.9,22 Following previous work,19 the PMFs were symmetrized by creating duplicate windows from the negative z domain and applying the WHAM method to all of the biased averaged ion density profiles (Figure 2B). Reference 19 shows that the expected errors in the barrier height computed with 2 ns sampling per window are roughly 0.3 kcal/mol. The error bars estimated here (0.1−0.2 kcal/mol) with the bootstrap method are consistent with those error estimates. The noted asymmetries suggest, however, that the PMF errors are somewhat larger due to longer time processes that are not completely sampled during the 4 ns trajectories. Thus, the PMF errors may be closer to 1 kcal/mol rather than the smaller errors suggested by the bootstrap method. The ion binding sites, located near the pore entrance on each side of the membrane, play an important role in the ion permeation process. The sites bind cations and reject anions.65 The previously reported outer binding site is at 11.3 Å, with an inner binding site at 9.7 Å.62 Our results are consistent with those previous results, and the binding sites are at approximately the same positions for the five temperatures examined (Figure 2B). The PMF for the case of a pure water bath possesses the same shape as for the water/glycerol bath. The binding site well depths and barrier height observed in the present study are slightly larger and smaller in magnitude, respectively, than those previously calculated.19 We note that the membrane lipid (here POPE and DMPC in ref 19) can have an effect on these quantities.33 Also, other studies9,19 explored the effect of the ion force field parameters on the PMFs. We do not pursue that issue here since we are mainly

the pure water bath case) are in decent agreement with previous simulations that employed the CHARMM27 force field (with the caveat, mentioned above, that our simulations involve a different lipid). Our results suggest that the free energy barrier in the water/ glycerol mixture is comparable to that for the pure water bath. The free energy barrier (water/glycerol bath) decreases from about 8.5 to 4.5 kcal/mol as the temperature is increased from 300 to 360 K. It can be noted that the PMF profile appears to be remarkably stable to increased temperatures up to 360 K, with a decrease in the barrier height at high temperatures (while retaining the same overall shape). The PMF at 345 K lies between the PMFs at 330 and 360 K, while the PMF at 315 K displays a barrier similar to that at 330 K; this result likely reflects uncertainties in the PMFs (discussed above). Below these PMFs will be scaled to smaller values (by a factor of 0.51) in order to obtain agreement with the experimental conductance estimate at 300 K. Decomposition into Enthalpy and Entropy Contributions. In order to gain further insight into the underlying mechanism of the decreasing free energy barrier height with increasing temperature, we computed separately the enthalpic and entropic parts of the free energy by numerical differentiation of the temperature-dependent free energy profiles (using the 300 and 360 K PMFs). The results are shown in Figure 3. (The presented results are computed directly from the raw PMF data in Figure 2 and are not scaled by the factor necessary to obtain agreement with experimental conductances.) It is clear that the free energy profile results from a significant compensation between the enthalpic and entropic contributions, since both Δhex and TΔsex are large and positive E



Article

though the protein environment is more constrained than the bulk solvent, structural changes occur at the pore center that lead to a positive USR value that well matches the value in bulk solution. (The change in the solvent reorganization term was also computed at z = ±5 Å, with a result of −4 kcal/mol; the result is still small in magnitude, but the negative value likely results from the more structured protein environment at these sites.) Also, using the qualitative Gaussian prediction of ⟨ε⟩/2 for the solvation free energy, these results suggest a barrier drop of 3.3 kcal/mol as the temperature increases from 300 to 360 K, consistent with the PMF result of 4 kcal/mol. It should be noted that a simple Gaussian model prediction for the barrier heights themselves (15.9, 14.2, and 12.7 kcal/mol) overestimates the barriers by roughly a factor of 2; this is due to neglected higher order contributions in the cumulant expansion that reflect interaction energy distribution asymmetries (see below). In addition, as discussed above, even if opposing errors of 1 kcal/mol occurred in the PMF barriers (leading to the 17.5 kcal/mol enthalpy barrier), the resulting prediction for the solvent reorganization energy change would be −10.5 kcal/mol, but the entropic portion of the barrier would still be positive with a magnitude of 11 kcal/mol. We next examine the fluctuation contributions. Figure 4 shows the log of the distributions of interaction energies of the

Figure 3. Enthalpy and entropy profiles of the K+ ion inside the gramicidin A channel at 330 K. The enthalpy in eq 1 and the entropy in eq 2 at 330 K were calculated from the PMFs at 300 and 360 K. The maximum size of the error bars for the enthalpic and entropic parts of the free energy is about 0.7 kcal/mol. The PMF, Δμ(cal) calculated by adding the enthalpic and entropic contributions agrees well with the PMF (Δμ(sim)) calculated from umbrella sampling simulations at 330 K (Figure 2). As discussed in the text, these enthalpy and entropy profiles were computed directly from the data of Figure 2 and were not scaled by the factor of 0.51 necessary for agreement with experimental conductance measurements.

at the pore center (28 and 22 kcal/mol, respectively). The enthalpic contribution is larger in magnitude, resulting in the free energy barrier. The relatively good agreement between the free energy profile obtained directly from the simulations and that constructed using eqs 1 and 2 suggests the numerical differentiation is relatively accurate. The estimated activation enthalpy (28 kcal/mol) is clearly too large in relation to experiment (the experimental estimate in ref 33 is 7.2 kcal/mol). When the result is scaled by the factor (0.51) necessary for agreement with the experimental conductance (below), a value of 14.3 kcal/mol results, still twice the experimental result (the lipid used in ref 33 is DiPhPC, and again this may also be a factor in the difference). First, if the PMF barrier height at 300 K were reduced by 1 kcal/mol, and the barrier at 360 K increased by 1 kcal/mol, the enthalpy barrier would be 17.5 kcal/mol, which leads to an estimate of 8.9 kcal/mol when scaled by 0.51, closer to experiment. This qualitative discussion shows that PMF errors on the order of 1 kcal/mol can show up as larger changes in the thermodynamic derivatives. The observation remains, however, that the PMF computed using the enthalpy and entropy values derived from data at 300 and 360 K matches well with the PMF from simulation (at 330 K), suggesting validity of the computed results. Second, the results may reflect limitations of the classical point charge force fields. It is more difficult to match experimental thermodynamic derivative quantities than the free energies themselves. Polarization effects may be important for obtaining results closer to experiment.25 To go further in the analysis, we employ the statistical mechanical theory discussed above. The computed interaction energy changes (Δ⟨ε⟩) for the ion moving from bulk solution to the pore center are 31.8, 28.4, and 25.3 kcal/mol at 300, 330, and 360 K, respectively. Taken together with the enthalpy change at 330 K computed from the temperature derivatives of the PMFs, this implies that the change in the solvent reorganization term (eq 9) is small in magnitude at the pore center (less than 0.5 kcal/mol). This indicates that, even

Figure 4. The log of the interaction energy distributions for the K+ ion located at z = −20 Å (bulk solvent) and z = 0 Å (center of the pore). The left three curves are for the ion in the bulk solvent at 300 K (black), 330 K (red), and 360 K (green). The right three curves are for the ion at the pore center at 300 K (blue), 330 K (cyan), and 360 K (brown). As discussed in the text, finite-size self-energy electrostatic corrections were not applied to the interaction energy distributions since they do not affect changes in mean interaction energy or fluctuation changes between the bulk solution and the pore.

ion with its surroundings in the bulk solvent and at the center of the pore. Several qualitative aspects of the distributions shed light on the ion interactions in the two environments. First, the distributions in the bulk water/glycerol solvent are relatively insensitive to temperature; this is reflected in consistent values both the mean binding energy and the fluctuation contribution to eq 7 (and resulting Gaussian free energy in the solvent) with increasing temperature. It can be noted that there is a slight positive skew in the distributions, evidenced by small positive third terms in the cumulant expansion (which are neglected in the Gaussian results). Second, the interaction energy distributions at the pore center are shifted substantially to higher energies, as discussed above. Third, the spreads of the F



Article

from the region approaching s = 0, but at values above the onset of the singularity, rather than in the range 10 < s < 30 (Figure 5A); this region produces correct diffusion constants in

distributions in the pore are significantly narrower than in the bulk solution; thus, the fluctuation contribution to the entropy change on entering the pore is relatively large and positive. Fourth, the distributions at the pore center exhibit a stronger temperature dependence than in the bulk solution, and especially at higher temperatures, they display a slight negative skew. Next, we focus on the Gaussian model results at 330 K for comparison with the entropy change computed directly from the PMFs. The numerical result for the entropy change from Figure 3 suggests a TΔsex value at the pore center of 22 kcal/ mol. Since the computed solvent reorganization change is small in magnitude, then the entropic contribution to eq 9 comes mainly from the fluctuation term. The Gaussian estimate for that term is 15.6 kcal/mol (with estimated errors on the order of 2 kcal/mol) or roughly 6.4 kcal/mol less than the value obtained from the PMFs. We suggest that the difference results from neglecting higher-order terms in the cumulant expansion, and the qualitative results of a positive skew in bulk solution and negative skew in the pore support that suggestion. (Equation 4 shows that the high-energy tail of the distribution yields the largest contribution to the free energy,55 so the positive skew in the bulk solution is the most important contributor to the deviation; the third-order term in the cumulant expansion in bulk solution is estimated to yield an entropic contribution of 6.7 kcal/mol, close to the above deviation.) The resulting Gaussian estimate for the free energy change is 12.8 kcal/mol, which when similarly shifted downward by 6.4 kcal/mol results in a value close to the observed 6 kcal/mol barrier obtained from the PMF. Thus, the Gaussian model leads to results consistent with those obtained from the PMFs. We note here that the interaction energy fluctuations were also examined at the locations z = ±5 away from the pore center, with an even larger reduction in fluctuations relative to the bulk, so it appears ion interaction energy fluctuations are reduced throughout the pore. We can draw several conclusions concerning the origin of the reduction of the PMF barrier height with increasing temperature. The drop of the barrier is reflected in the positive entropy change for the ion entering the pore (from eq 2). The above analysis indicates that the positive entropy change in turn is predominantly due to decreased fluctuations in the ion interaction energy upon entering the pore, with a smaller contribution from the change in the solvent reorganization energy. These results suggest a delicate balance of protein structure and fluctuation; the protein rearranges due to the presence of the ion in a way that balances the solvent reorganization energy in the bulk solution, yet the structure of the protein limits the ion interaction energy fluctuations. The general consistency of the Gaussian model results with those obtained from the computed PMFs lends independent support to the validity of the PMF-derived enthalpy and entropy profiles. Diffusion Constants. During umbrella sampling for the PMF calculations, a biased systematic force was applied to the K+ ion, so diffusion coefficients calculated directly from the MSD and VACF yield incorrect results. One-dimensional diffusion constants for the K+ ion in the gramicidin channel were thus calculated using the GLE-HO method. Consistent with the results in ref 60, we observe that D̂ (s) diverges near s = 0 due to numerical artifacts resulting from finite sampling. The location of this singularity depends on the sampling time.60 Diffusion coefficients were calculated by linearly extrapolating

Figure 5. Diffusion coefficient profile for the K+ ion inside the HD form at 360 K. (A) The y-intercept of the red line represents the diffusion coefficient of K+ in the bulk region, calculated by linear fitting from the region 1 < s < 3. (B) Relative diffusion coefficients were calculated from 81 windows and symmetrized. The diffusion coefficient in the bulk region was set to 1. The red line represents sigmoidal fitting.

comparison with the bulk ion data given above. An alternative theoretical approach involving position correlation functions and Bayesian analysis has been developed by Hummer;66 that method was found to produce results consistent with the approach used here. Finally, we calculated the z-dependent relative diffusion constant in the channel in relation to the bulk value (Figure 5B). The diffusion constant of the ion in the bulk region for the HD form increases as temperature increases: 0.09 (Å2/ps) at 300 K, 0.14 (Å2/ps) at 330 K, and 0.22 (Å2/ps) at 360 K. There are some fluctuations of diffusion constants inside the gramicidin channel, so relative diffusion coefficients are symmetrized by averaging both sides, followed by sigmoidal fitting (Figure 5B). Increasing temperature did not dramatically change the ratio of the fitted diffusion coefficient inside the channel to the bulk region: 28% at 300 K, 23% at 330 K, and 19% at 360 K. Our results are similar to the results from other groups,19 although the diffusion coefficients are somewhat smaller inside the pore due to fitting at smaller s values, as suggested in ref 60. Edwards et al.67 showed that a diffusion coefficient inside the channel that is 20% of the bulk value fits the linear I−V curve obtained from Brownian dynamics simulations and matches experimental data. We chose the Langevin damping coefficient as 5 ps−1 for the NPT simulations. It should be noted that the choice of Langevin G



Article

HOLE method.70 The sampling distance between planes was taken as 0.25 Å, and a termination radius value of 5 Å was used. Position data for the pore radius calculations were output every 10 ps and averaged. The lines with filled circles are for the radius calculated by holding an ion at a fixed location during umbrella sampling (Figure 6). The lines (without symbols) are

damping coefficient can in principle affect the magnitude of the diffusion coefficient. When we changed its magnitude from 5 to 1, however, the ratio of diffusion coefficient of K+ ion in the bulk to the center of the channel changed by less than 10%. The Maximum Conductance. The anticipated potassium ion conductances from the experiments of Cuppoletti et al.28−30 (obtained by scaling the measured proton conductance downward by a factor of 30) are 4.4 pS at 300 K and 8.8 pS at 360 K. As mentioned above, these K+ conductance estimates are smaller than the K+ data presented in refs 33 and 35 (roughly 24.5 pS at 300 K), and this is likely due to the different membrane phospholipids and solvent mixtures employed in the present study. The experiments in ref 33 were performed with a 1 M concentration of KCl, and we estimate the maximum conductance at saturation to be 6.6 pS (see ref 34, p 367). Using the raw PMF data, the calculated maximum conductance of the HD form increases from 0.007 to 1.19 pS as the temperature increases from 300 to 360 K. In terms of the conductance calculation, the strong increase is due to (1) the decrease in the barrier height, (2) the scaling of the PMF in the exponential of eq 16 by the temperature, and (3) the increased diffusion constants. Thus, the conductances obtained by using the computed PMFs obtained from the simulations are not in agreement with experiment, as has been noted before.19 In that study, corrections to the PMFs were computed that may arise from finite system sizes and dielectric effects in the membrane, which reduced the PMF barriers substantially. Rather than computing similar corrections to the PMFs, we alternatively choose to scale the PMF at 300 K in order to obtain agreement with the experimental conductance. We found that scaling the PMF by a factor of 0.51 resulted in good agreement. The PMF barrier height with the scaling is then 4.3 kcal/mol. All PMF curves at the other temperatures were similarly scaled by the same factor. Then the barriers at 330 and 360 K decreased to 3.2 and 2.2 kcal/mol, respectively. The resulting computed conductances are 6.6 pS at 300 K, 31.9 pS at 330 K, and 67.3 pS at 360 K, showing substantial deviation from experiment at high temperature (the estimated experimental result at 360 K is 13.2 pS). The ratio of the computed high and low temperature conductances is 10, significantly larger than the observed factor of 2 for proton conduction. We note that, if the K+ data in Figure 5 of ref 33 are extrapolated to 360 K, the predicted conductance ratio is roughly 7, in decent agreement with our calculated estimate. The deviation of our calculated conductance ratio for K+ transport from the experimental ratio for proton conduction may be linked both to short-time blockages of the pore at higher temperatures and to fundamental differences between proton and potassium ion transport through the pore (below). We also computed the conductance at 360 K using the PMF profile at 300 K. The result is 15.1 pS, yielding a conductance ratio of 2.3, well below the value of 7 obtained above by extrapolating the K+ data of ref 33. We suggest that this result lends independent support to the notion that the free energy barrier drops as the temperature increases, in agreement with the previous discussion of Urry et al.35 This can be taken as a further indication that the entropy change for the ion entering the pore is positive. Pore Radius. The pore radius can play a significant role in ion permeation since it may affect the motion of the ion and the pore waters.68,69 In order to investigate the structural changes in the gramicidin pore, we calculated the pore radius using the

Figure 6. Pore radius profiles for the gramicidin A channel. The lines with filled circles indicate the radius was calculated by holding an ion at a fixed point during umbrella sampling. The lines without symbols are for the case of no ions in the pore during 3 ns simulations.

for the pore radius calculated without an ion inside the pore during a 3 ns simulation. The ionic radius of K+ is 1.37 Å.71 The radius of the ion-filled pore is smaller than the radius of the ion except at the channel center (dimer interface). The ion-free pore radius profile clearly exhibits some blockage near the channel entrances and a larger free volume along the pore and near the channel center at the lower temperatures. Interestingly, for the ion-free form we observed a drop in the average pore radius at the center of the channel at 360 K. The formyl group at the N-terminus frequently blocked the center of the channel. Typically the channel is filled with around 8 water molecules, and those water molecules are hydrogen bonded to each other. But when the center of the HD structure is occupied by a formyl group, the water chain is broken. Conversely, when the formyl group receded from the dimer interface, the water chain was again restored. This observation may help to reconcile the calculated ratio of 10 for the conductances at 360 and 300 K with the experimental ratio of about 2. Frequent short-term blockages of the unoccupied pore by the formyl groups may lead to lower average conductances at the high temperature due to the breaking of the water wire in the pore. The process by which we calculate the conductance from the PMF (with an ion constrained to small windows along the z direction) does not take into account the possibility of blocks of the unoccupied channel by the formyl groups since the ions are inserted into the pore. Therefore, the conductances computed from the PMFs and diffusion theory are expected to result in overestimates. The observed time scale for the formyl blocks is on the order of 100 ps, which would not show up in single-channel conduction traces; rather, the measured openstate conductance is likely an average that includes the formyl blockages. The pore radius when the ion is held at a restrained position is significantly less than the pore radius calculated without an ion inside the channel and is largely independent of temperature. These results illustrate the impact of the ion on H



Article

and 5. The average distance is 2.15 Å at 300 K, 2.19 Å at 330 K, and 2.43 Å at 360 K. Next, we defined an axis for each monomer as the line which connects the center of mass of the backbone oxygens of residues 1, 3, 5 and 8, 10, 12. The average bending angle increases from 4.7° to 7.0° as the temperature increases from 300 to 360 K. Again, the results indicate that high temperatures lead to somewhat larger fluctuations of the channel.

the local peptide environment relative to the ion-free case. The main reason for the decrease in the pore radius is the strong attractive force between the ion and the gramicidin backbone carbonyls.72 The pore radius at the center displays a maximum because there are no backbone oxygens at that position. Our results agree with those previously reported by Siu and Böckmann.9 There appears to be little temperature dependence to the pore radius around sites occupied by the ions, unlike for the unoccupied HD case. The Number of Hydrogen Bonds at the Junction. It is known that the mechanism of gating of the HD form is related to association and dissociation of the dimer.34,73 Recently, Lu suggested that the HD form can be in a closed intermediate state without a large separation of the dimer.74 The number of hydrogen bonds between dimers likely plays an important role in opening and closing of the channel. We investigated the distribution of the number of hydrogen bonds at the interface at 300, 330, and 360 K using the VMD code. The following conditions for the existence of H-bonds were used: 3.2 Å for donor−acceptor distance and 30° for the angle cutoff. There is a maximum of six hydrogen bonds at the junction (H1−O5, O1−H5, H3−O3, O3−H3, H5−O1, O5−H1; see ref 73 for definition of the H-bonds). During a 21 ns simulation, the number of H-bonds fluctuated widely; the distribution of Hbonds is shown in Figure 7. As the temperature increases, the

■

CONCLUSIONS The present modeling study was motivated by the recent measurements of proton conduction through the gramicidin channel at temperatures approaching the boiling point of water28−30 and other studies of the temperature-dependent transport of K+ 33,35 and H+.33 In this initial study, we chose to model K+ transport due to the severe difficulty in accurately representing the chemical aspects of proton migration. The experiments in refs 33 and 35 did not report results at the highest temperature examined here, so we scaled the measured proton conductance results based on the data for both ions in ref 33 at 300 K. The mechanisms of K+ and H+ transport surely differ, and thus the temperature dependences may also differ. Both ions display activated transport, however, with the enthalpy barrier for K+ (7.2 kcal/mol) being twice that for H+ (3.7 kcal/mol).33 The main conclusions of our paper are first that, in agreement with experiment, the channel remains stable over 100 ns time scales at 360 K and second that the free energy barrier height for K+ transport drops with increasing temperature, implying a positive entropy change for the ion entering the pore. There is a large compensation between the enthalpy and entropy parts of the free energy, leading to a free energy barrier with a significantly smaller magnitude than the contributing parts. The computed conductance ratios at 360 and 300 K, along with data from ref 33 extrapolated to 360 K, are in the range 7−10. When the computed PMF at 300 K is used for a conductance calculation at 360 K, however, a much smaller ratio is obtained. This provides further, although indirect, evidence of a drop in the free energy barrier with increasing temperature. By employing a statistical mechanical theory for the enthalpy and entropy, we found that the solvent reorganization energy change upon entering the channel pore is small and that the positive entropy change is due mainly to reduced fluctuations in the ion interaction energy in the pore relative to those in the bulk solution. The reduced fluctuations are in turn due to the more structured environment in the pore, although there are still clearly substantial pore fluctuations that allow for rapid ion transit through the channel. (The computational studies of Bastug et al.75 and Roux and Karplus76 show the importance of the protein flexibility in ion transport.) The observation that the positive entropy change is due to reduced fluctuations may at first appear counterintuitive. The ion entropy profile is a partial molar quantity, however, and what the results suggest is that the ion causes a greater restriction on the nearby solvent molecules when inserted into the bulk solution relative to the change in the surroundings upon inserting the ion in the pore (which is already somewhat restricted by the protein structure). The values computed for the enthalpy and entropy change are too large in magnitude, even when the PMF barriers are scaled downward by a factor chosen to yield agreement with experimental conductances. The downward scaling presumably mimics effects of polarization in the peptide−membrane

Figure 7. Probability distributions of H-bonds at the channel junction during a 21 ns simulations. The data were divided into seven blocks, where one block corresponds to 3 ns simulations, and the probability distributions and error bars were calculated from the seven blocks.

distribution shifts to smaller numbers of H-bonds. No complete dissociations of the dimer were observed in our simulations, however. These results indicate the larger fluctuations of the channel at high temperatures. We also note that these simulations imply that the overall structure of the channel is stable over rather long time scales at temperatures up to 360 K. In order to further test the long-time stability of the channel, one simulation was run for 100 ns at 360 K, and the dimeric structure was again preserved. Distance and Bending Angle Difference between Monomers. In order to further investigate the effect of temperature on the channel structure, we computed the average distance between monomers and the average bending angle difference between monomers during a 3 ns simulation. The distance between the monomers is defined as the distance of the center of mass for the backbone oxygens for residues 1, 3, I



Article

It is hoped that the results discussed in this paper will stimulate further experimental and theoretical work aimed at a more thorough understanding of the temperature dependence of ion transport through a range of ion channels. Further, the results are consistent with the experimental observation that the gramicidin A channel can function at surprisingly high temperatures and may thus find use in device applications such as fuel cell technology. A more speculative suggestion is that, by controlling the rigidity of synthetic channels, the thermodynamics of transport may be tunable for particular device applications at elevated temperatures.

environment. The deviation of the enthalpy barrier from experiment may be tied to limitations of the classical force fields; thermodynamic derivative quantities are more difficult to accurately model than are the free energies themselves. Related classical models have been shown to yield bulk hydration entropies in relatively good agreement with experiment,58 however, so this suggests that if there are force field deficiencies, they are most important for the ion−protein interactions. Nevertheless, we suggest that the observed positive entropy change likely is a common feature for ions entering narrow pores in proteins. The positive entropy change observed in the present study of K+ transport differs from the experimental findings for H+ transport in ref 33, where negative entropies were reported. First, as discussed above, the underlying mechanisms of K+ and H+ transport differ and thus the temperature dependences may differ. Second, the analysis of ref 33 relied on Eyring rate theory; their results for the enthalpy barrier are unambiguous, but the computed entropy change involves the prefactor in the rate expression. It has been suggested previously77 that transition state theory is not necessarily applicable to ion transport through channels with barriers of modest size. The present results also differ from a previous computational study of K+ ion enthalpy and entropy profiles through model channels.53 In that study, restraints were placed on the protein in order to maintain a well-defined structure. Positive enthalpies and negative entropies were reported for ion entry into the pore. On the basis of the above discussion, we suggest that placement of the constraints on the protein results in a significantly negative ΔUSR (solvent reorganization energy change) value. The separate USR values in the bulk solvent and protein are positive and large in magnitude, and placement of the constraints will likely reduce the USR value when the ion is in the pore, resulting in the observed negative entropy change (and less positive enthalpy change). Of course, the constraints could also alter the fluctuation part of the entropy, but the results presented here suggest the change is of smaller magnitude than that for the solvent reorganization term. Finally, frequent blocks of the unoccupied pore by the terminal formyl groups near the dimer interface were observed at 360 K. This observation may help to rationalize the small experimentally measured factor of 2 increase of proton conductance between 300 and 360 K. Extrapolation of the conductance data from ref 33 to higher temperatures (and our computed conductances) implies a significantly larger increase in the conductance with temperature; the formyl blockage may effectively reduce the measured currents at the higher temperatures. The time scale for the formyl blocks is much shorter than what could be observed in patch clamp experiments. In agreement with an early pioneering study of the gramicidin A channel,76 the present study suggests that entropic effects are important for single-ion motion through narrow channels. It is interesting that previous computational studies have suggested that entropic effects are relatively small for the K+/Na+ filter contribution to the selectivity in the potassium channel;78−81 those results may indicate a delicate cancellation of nonzero entropy change values for the two different ions moving from bulk solution into the pore, although we emphasize that the entropy change in bulk solution also contributes to the overall selectivity (the K+/Na+ Tsex change in bulk water is −2.6 kcal/mol82).

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address §

H.D.S.: Department of Physics, Indiana University−Purdue University, Indianapolis, Indianapolis, IN 46202. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank John Cuppoletti, Frank Pinski, Manori Jayasinghe, and Dilip Asthagiri for helpful discussions. We gratefully acknowledge a generous grant of computer time at the Ohio Supercomputer Center and financial support from the National Science Foundation (CHE-0709560 and CHE-1011746) and the University of Cincinnati Research Council.

■

REFERENCES

(1) Burkhart, B. M.; Li, N.; Langs, D. A.; Pangborn, W. A.; Duax, W. L. The Conducting Form of Gramicidin A is a Right-Handed DoubleStranded Double Helix. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 12950− 12955. (2) European Pharmacopoeia; Council of Europe: Strasbourg, France, 2004. (3) Urry, D. W. The Gramicidin A Transmembrane Channel: A Proposed π(L,D) Helix. Proc. Natl. Acad. Sci. U. S. A. 1971, 68, 672−676. (4) Veatch, W. R.; Fossel, E. T.; Blout, E. R. The Conformation of Gramicidin A. Biochemistry 1974, 13, 5249−5256. (5) Glowka, M. L.; Olczak, A.; Bojarska, J.; Szczesio, M.; Duax, W. L.; Burkhart, B. M.; Pangborn, W. A.; Langs, D. A.; Wawrzak, Z. Structure of Gramicidin D-RbCl Complex at Atomic Resolution from LowTemperature Synchrotron Data: Interactions of Double-Stranded Gramicidin Channel Contents and Cations with Channel Wall. Acta Crystallogr. 2005, D61, 433−441. (6) Wallace, B. A. Gramicidin Channels and Pores. Annu. Rev. Biophys. Biophys. Chem. 1990, 19, 127−157. (7) Killian, J. A.; Prasad, K. U.; Hains, D.; Urry, D. W. The Membrane as an Environment of Minimal Interconversion. A Circular Dichroism Study on the Solvent Dependence of the Conformational Behavior of Gramicidin in Diacylphosphatidylcholine Model Membranes. Biochemistry 1988, 27, 4848−4855. (8) Lin, T.-H.; Huang, H.-b.; Wei, H.-A.; Shiao, S.-H.; Chen, Y. C. The Effect of Temperature and Lipid on the Conformational Transition of Gramicidin A in Lipid Vesicles. Biopolymers 2005, 78, 179−186. (9) Siu, S. W. I.; Böckmann, R. A. Low Free Energy Barrier for Ion Permeation Through Double-Helical Gramicidin. J. Phys. Chem. B 2009, 113, 3195−3202. (10) Rawat, S. S.; Kelkar, D. A.; Chattopadhyay, A. Monitoring Gramicidin Conformations in Membranes: A Fluorescence Approach. Biophys. J. 2004, 87, 831−843.

J



Article

(11) Kelkar, D. A.; Chattopadhyay, A. The Gramicidin Ion Channel: A Model Membrane Protein. Biochim. Biophys. Acta, Biomembr. 2007, 1768, 2011−2025. (12) Olah, G. A.; Huang, H. W.; Liu, W. H.; Wu, Y. L. Location of Ion-Binding Sites in the Gramicidin Channel by X-Ray Diffraction. J. Mol. Biol. 1991, 218, 847−858. (13) Andersen, O. S.; Apell, H. J.; Bamberg, E.; Busath, D. D.; Koeppe, R. E., II; Sigworth, F. J.; Szabo, G.; Urry, D. W.; Woolley, A. Gramicidin Channel Controversy−the Structure in a Lipid Environment. Nat. Struct. Biol. 1999, 6, 609. (14) Cross, T. A.; Arseniev, A.; Cornell, B. A.; Davis, J. H.; Killian, J. A.; Koeppe, R. E., II; Nicholson, L. K.; Separovic, F.; Wallace, B. A. Gramicidin Channel Controversy−Revisited. Nat. Struct. Biol. 1999, 6, 610−611. (15) Burkhart, B. M.; Duax, W. L. Gramicidin Channel Controversy Reply. Nat. Struct. Biol. 1999, 6, 611−612. (16) Dzikovski, B. G.; Borbat, P. P.; Freed, J. H. Channel and Nonchannel Forms of Spin-Labeled Gramicidin in Membranes and Their Equilibria. J. Phys. Chem. B 2011, 115, 176−185. (17) Allen, T. W.; Bastug, T.; Kuyucak, S.; Chung, S. H.; Gramicidin, A. Channel as a Test Ground for Molecular Dynamics Force Fields. Biophys. J. 2003, 84, 2159−2168. (18) Allen, T. W.; Andersen, O. S.; Roux, B. Energetics of Ion Conduction Through the Gramicidin Channel. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 117−122. (19) Allen, T. W.; Andersen, O. S.; Roux, B. Ion Permeation Through a Narrow Channel: Using Gramicidin to Ascertain All-Atom Molecular Dynamics Potential of Mean Force Methodology and Biomolecular Force Fields. Biophys. J. 2006, 90, 3447−3468. (20) Ingólfsson, H. I.; Li, Y.; Vostrikov, V. V.; Gu, H.; Hinton, J. F.; Koeppe, R. E., II; Roux, B.; Andersen, O. S.; Gramicidin, A. Backbone and Side Chain Dynamics Evaluated by Molecular Dynamics Simulations and Nuclear Magnetic Resonance Experiments. I: Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 7417− 7426. (21) Bastug, T.; Kuyucak, S. Test of Molecular Dynamics Force Fields in Gramicidin A. Eur. Biophys. J. 2005, 34, 377−382. (22) Bastug, T.; Patra, S. M.; Kuyucak, S. Molecular Dynamics Simulations of Gramicidin A in a Lipid Bilayer: From StructureFunction Relations to Force Fields. Chem. Phys. Lipids 2006, 141, 197−204. (23) Mustafa, M.; Busath, D. D. The Gramicidin Channel Ion Permeation Free-Energy Profile: Direct and Indirect Effects of CHARMM Force Field Improvements. Interdiscip. Sci. 2009, 1, 113−127. (24) De Fabritiis, G.; Coveney, P. V.; Villa-Freixa, J. Energetics of K+ Permeability through Gramicidin A by Forward-Reverse Steered Molecular Dynamics. Proteins: Struct. Funct. Bioinf. 2008, 73, 185−194. (25) Patel, S.; Davis, J. E.; Bauer, B. A. Exploring Ion Permeation Energetics in Gramicidin A Using Polarizable Charge Equilibration Force Fields. J. Am. Chem. Soc. 2009, 131, 13890−13891. (26) Bastug, T.; Patra, S. M.; Kuyucak, S. Finite System and Periodicity Effects in Free Energy Simulations of Membrane Proteins. Chem. Phys. Lett. 2006, 425, 320−323. (27) Myers, V. B.; Haydon, D. A. Ion Transfer across Lipid Membranes in the Presence of Gramicidin A: II. The Ion Selectivity. Biochim. Biophys. Acta, Biomembr. 1972, 274, 313−322. (28) Korfhagen, S. Stabilization of Scaffold-Supported, Photopolymerized Bilayer Lipid Membranes with Gramicidin-D for Novel Fuel Cells. M.Sc. Thesis, University of Cincinnati, 2008. (29) Cuppoletti, J.; Birn, S. S.; Tewari, K. P.; Chakrabarti, J.; Malinowska, D. H. Studies of Ion Transport at Extreme Temperature. FASEB J. 2008, 22, 1201.24. (30) Cuppoletti, J.; Malinowska, D. H. In Advances in Composite Materials − Analysis of Natural and Man-Made Materials; Tesinova, P., Ed.; Intech: Rijeka, 2011; Chapter 10, pp 285−294. (31) Smitha, B.; Sridhar, S.; Khan, A. A. Solid Polymer Electrolyte Membranes for Fuel Cell Applications − a Review. J. Membr. Sci. 2005, 259, 10−26.

(32) Qin, Z.; Tepper, H. L.; Voth, G. A. Effect of Membrane Environment on Proton Permeation through Gramicidin A Channels. J. Phys. Chem. B 2007, 111, 9931−9939. (33) Chernyshev, A.; Cukierman, S. Thermodynamic View of Activation Energies of Proton Transfer in Various Gramicidin A Channels. Biophys. J. 2002, 82, 182−192. (34) Hille, B. Ion Channels of Excitable Membranes; Sinauer Associates Inc.: Sunderland, MA, 2001. (35) Urry, D. W.; Alonso-Romanowski, S.; Venkatachalam, C. M.; Bradley, R. J.; Harris, R. D. Temperature Dependence of Single Channel Currents and the Peptide Libration Mechanism for Ion Transport through the Gramicidin A Transmembrane Channel. J. Membr. Biol. 1984, 81, 205−217. (36) MacKerell, A. D.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (37) Townsley, L. E.; Tucker, W. A.; Sham, S.; Hinton, J. F. Structures of Gramicidins A, B, and C Incorporated into Sodium Dodecyl Sulfate Micelles. Biochemistry 2001, 40, 11676−11686. (38) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (39) Mustafa, M.; Busath, D. D. The Gramicidin Channel Ion Permeation Free-Energy Profile: Direct and Indirect Effects of CHARMM Force Field Improvements. Interdiscip. Sci. Comput. Life Sci. 2009, 1, 113−127. (40) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613−4621. (41) Siu, S. W.; Böckmann, R. A. Electric Field Effects on Membranes: Gramicidin A as a Test Ground. J. Struct. Biol. 2007, 157, 545−556. (42) Rand, R. P.; Fuller, N.; Parsegian, V. A.; Rau, D. C. Variation in Hydration Forces between Neutral Phospholipid Bilayers: Evidence for Hydration Attraction. Biochemistry 1988, 27, 7711−7722. (43) Tieleman, D. P.; Berendsen, H. Molecular Dynamics Simulations of a Fully Hydrated Dipalmitoylphosphatidylcholine Bilayer with Different Macroscopic Boundary Conditions and Parameters. J. Chem. Phys. 1996, 105, 4871−4880. (44) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N·log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (45) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (46) Kalé, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283−312. (47) Walton, E. B.; Vanvliet, K. J. Equilibration of Experimentally Determined Protein Structures for Molecular Dynamics Simulation. Phys. Rev. E 2006, 74, 061901. (48) Chipot, C.; Pohorille, A. Free Energy Calculations; Springer: New York, 2007. (49) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I: The Method. J. Comput. Chem. 1992, 13, 1011−1021. (50) Roux, B. The Calculation of the Potential of Mean Force Using Computer Simulations. Comput. Phys. Commun. 1995, 91, 275−282. (51) Grossfield, A. WHAM code: http://membrane.urmc.rochester. edu/content/wham/. (52) Hinsen, K.; Roux, B. Potential of Mean Force and Reaction Rates for Proton Transfer in Acetylacetone. J. Chem. Phys. 1997, 106, 3567−3577. (53) Portella, G.; Hub, J. S.; Vesper, M. D.; de Groot, B. L. Not Only Enthalpy: Large Entropy Contribution to Ion Permeation Barriers in Single-File Channels. Biophys. J. 2008, 95, 2275−2282. K



Article

(54) Beck, T. L.; Paulaitis, M. E.; Pratt, L. R. The Potential Distribution Theorem and Models of Molecular Solutions; Cambridge University Press: Cambridge, 2006. (55) Beck, T. L. Hydration Free Energies by Energetic Partitioning of the Potential Distribution Theorem. J. Stat. Phys. 2011, 145, 335−354. (56) Asthagiri, D.; Pratt, L. R.; Paulaitis, M. E. Role of Fluctuations in a Snug-Fit Mechanism of KcsA Channel Selectivity. J. Chem. Phys. 2006, 125, 024701. (57) Ben-Amotz, D.; Underwood, R. Unraveling Water’s Entropic Mysteries: A Unified View of Nonpolar, Polar, and Ionic Hydration. Acc. Chem. Res. 2008, 41, 957−967. (58) Beck, T. L. A Local Entropic Signature of Specific Ion Hydration. J. Phys. Chem. B 2011, 115, 9776−9781. (59) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (60) Mamonov, A. B.; Kurnikova, M. G.; Coalson, R. D. Diffusion Constant of K+ inside Gramicidin A: A Comparative Study of Four Computational Methods. Biophys. Chem. 2006, 124, 268−278. (61) Levitt, D. G. Interpretation of Biological Ion Channel Flux Data−Reaction-Rate versus Continuum Theory. Annu. Rev. Biophys. Biophys. Chem. 1986, 15, 29−57. (62) Allen, T. W.; Andersen, O. S.; Roux, B. Molecular Dynamics − Potential of Mean Force Calculations as a Tool for Understanding Ion Permeation and Selectivity in Narrow Channels. Biophys. Chem. 2006, 124, 251−267. (63) Roux, B. Statistical Mechanical Equilibrium Theory of Selective Ion Channels. Biophys. J. 1999, 77, 139−153. (64) Hollerbach, U.; Chen, D.-P.; Eisenberg, R. S. Two- and ThreeDimensional Poisson-Nernst-Planck Simulations of Current Flow Through Gramicidin A. J. Sci. Comput. 2001, 16, 373−409. (65) Bastug, T.; Kuyucak, S. Energetics of Ion Permeation, Rejection, Binding, and Block in Gramicidin A from Free Energy Simulations. Biophys. J. 2006, 90, 3941−3950. (66) Hummer, G. Position-Dependent Diffusion Coefficients and Free Energies from Bayesian Analysis of Equilibrium and Replica Molecular Dynamics Simulations. New J. Phys. 2005, 7, 34. (67) Edwards, S.; Corry, B.; Kuyucak, S.; Chung, S. H. Continuum Electrostatics Fails to Describe Ion Permeation in the Gramicidin Channel. Biophys. J. 2002, 83, 1348−1360. (68) Tang, Y. Z.; Chen, W. Z.; Wang, C. X. Molecular Dynamics Simulations of the Gramicidin A-Dimyristoylphosphatidylcholine System with an Ion in the Channel Pore Region. Eur. Biophys. J. 2000, 29, 523−534. (69) Chiu, S. W.; Subramaniam, S.; Jakobsson, E.; McCammon, J. A. Water and Polypeptide Conformations in the Gramicidin Channel. A Molecular Dynamics Study. Biophys. J. 1989, 56, 253−261. (70) Smart, O. S.; Goodfellow, J. M.; Wallace, B. A. The Pore Dimensions of Gramicidin A. Biophys. J. 1993, 65, 2455−2460. (71) Shannon, R. D. Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides. Acta Crystallogr. 1976, A32, 751−767. (72) Tang, Y. Z.; Chen, W. Z.; Wang, C. X. Molecular Dynamics Simulations of the Gramicidin A-Dimyristoylphosphatidylcholine System with an Ion in the Channel Pore Region. Eur. Biophys. J. 2000, 29, 523−534. (73) Miloshevsky, G. V.; Jordan, P. C. Gating Gramicidin Channels in Lipid Bilayers: Reaction Coordinates and the Mechanism of Dissociation. Biophys. J. 2004, 86, 92−104. (74) Lu, H. P. Probing Single-Molecule Protein Conformational Dynamics. Acc. Chem. Res. 2005, 38, 557−565. (75) Bastug, T.; Gray-Weale, A.; Patra, S. M.; Kuyucak, S. Role of Protein Flexibility in Ion Permeation: A Case Study in Gramicidin A. Biophys. J. 2006, 90, 2285−2296. (76) Roux, B.; Karplus, M. Ion Transport in a Model Gramicidin Channel: Structure and Thermodynamics. Biophys. J. 1991, 59, 961− 981. (77) Roux, B.; Karplus, M. Ion Transport in a Gramicidin-Like Channel: Dynamics and Mobility. J. Phys. Chem. 1991, 95, 4856−4868.

(78) Roux, B.; Bernache, S.; Egwolf, B.; Lev, B.; Noskov, S. Y.; Rowley, C. N.; Yu, H. Ion Selectivity in Channels and Transporters. J. Gen. Physiol. 2011, 137, 415−426. (79) Noskov, S. Y.; Roux, B. Ion Selectivity in Potassium Channels. Biophys. Chem. 2006, 124, 279−291. (80) Dixit, P. D.; Merchant, S.; Asthagiri, D. Ion Selectivity in the KcsA Potassium Channel from the Perspective of the Ion Binding Site. Biophys. J. 2009, 96, 2138−2145. (81) Varma, S.; Rogers, D. M.; Pratt, L. R.; Rempe, S. B. Design Principles for K+ Selectivity in Membrane Transport. J. Gen. Physiol. 2011, 137, 479−488. (82) Schmid, R.; Miah, A. M.; Sapunov, V. N. A New Table of the Thermodynamic Quantities of Ionic Hydration: Values and Some Applications (Enthalpy−Entropy Compensation and Born Radii). Phys. Chem. Chem. Phys. 2000, 2, 97−102.

L


Temperature Dependence of Gramicidin Channel Transport and

Recommend Documents