5318
J. Phys. Chem. B 2009, 113, 5318–5326
Hydrogen Bonds in Membrane Proteins Sheh-Yi Sheu,*,†,‡ Edward W. Schlag,*,‡,§ Heinrich L. Selzle,§ and Dah-Yen Yang*,‡,| Department of Life Sciences and Institute of Biomedical Informatics, National Yang-Ming UniVersity, Taipei 112, Taiwan, Institut fu¨r Physikalische and Theoretische Chemie, TU-Mu¨nchen, Lichtenbergstrasse 4, 85748 Garching, Germany, and Institute of Atomic and Molecular Science, Academia Sinica, Taipei 106, Taiwan ReceiVed: December 8, 2008; ReVised Manuscript ReceiVed: February 16, 2009
Hydrogen bonds are essential tie points inside protein structures. They undergo dynamic rupture and rebonding processes on the time scale of tens of picoseconds. Proteins can partially rearrange during such ruptures. In previous work, we performed molecular dynamics simulations of these fluctuating hydrogen bonds. This indicated long-range entropy and energy contributions extending far into the liquid environment. The results showed that the binding of a given hydrogen bond is much reduced as a result of these interactions in water, as is required for biological activity and in very good confirmation of known experimental results. The larger water environment directly interacts with the hydrogen bond essentially due to long-range molecular interactions. Such a substantial lowering of the energy of the hydrogen bond in water brings it into the range of activation by many biological processes (Sheu et al. Chem. Phys. Lett. 2008, 462, 1-5). Thus, the water medium profoundly increases the rate. Furthermore, very large entropic changes are associated with the rupture of hydrogen bonds in water, whereas no such effects are seen for the isolated molecule. Interestingly, such an increase in rates in water is still accompanied by a large negative change in entropy in the extended solvent environment, and this reduces the rate by some 2 orders of magnitude. Recent molecular dynamics experiments in D2O substantiate this model and show a large solvent isotope effect. In this work, we used lipids as the environment for the hydrogen bond and discovered that the energy is also reduced from that found in the isolated molecule, but not as far as in water. On the other hand, we found that no entropy penalty exists for breaking the hydrogen bond in lipids, as seen for water. These two effects compensate, even though the energy is some 2 times larger. The entropic penalty is reduced such that the rate is higher than in water despite the higher energy. This is a significant result for understanding the reactivity and dynamics of proteins in lipids. It should be noted that these are very important solvent effects on entropies and free energies that are not usually reflected in statistical thermodynamic computations for reactants and products. The very longrange effect of the solvent makes substantial contributions to kinetic rate constants and is readily evaluated in this kinetic method. To ignore these long-range environmental effects on the entropy can lead to very spurious results when calculating rates of protein mobilities. Hence, the results not only agree very well with the known hydrogen-bond energies directly as a result of various environmental factors, but even correctly predict a phase transition in the lipid. I. Introduction 1-6
The hydrogen bond is one of the key elements determining the dynamics and stability of proteins.7,8 Its strength and directionality affect the property of protein mobility in solvated systems.9-12 Furthermore, such hydrogen bonds might not be fixed, but might instead be made and broken on a time scale tens of picoseconds. This permits temporary tie points to be formed during the folding process. An isolated hydrogen bond is some 6 kcal in strength;13-15 this is lowered to some 1.6 kcal in the water medium.16-21 We showed in previous work that this strong change is mostly due to the direct solvation of the hydrogen bond by the water medium, with concomitant large changes in the long-range water structure upon activation.22 Even though the rate is enhanced by the low activation energy in * To whom correspondence should be addressed. E-mail: sysheu@ ym.edu.tw (S.-Y.S.),
[email protected] (E.W.S.),
[email protected]. edu.tw (D.-Y.Y.). † National Yang-Ming University. ‡ S.-Y.S., E.W.S., and D.-Y.Y. contributed equally to this work. § TU-Mu¨nchen. | Academia Sinica.
water, the rate is inhibited by 2 orders of magnitude because of a strong entropic penalty upon activation. Such entropic effects are not contained in the usual partition functions of the molecular system but come from bulk effects of the medium, and as such, they assume an additional dominant importance for all rate processes. Hence, the hydrogen bond should be understood here not as a molecular property, but as a holistic property of the fluctuating medium interacting with a weak hydrogen bond. Any change in the medium strongly affects the entropic bottleneck, despite its low energy.17 In this work, we investigated the effects of changing the environment of a membrane protein that is embedded in a heterogeneous environment consisting of the phospholipid bilayer membrane regions with water surrounding the extramembrane domains. Typically, the water-solvated part of the membrane protein has a diverse structure, whereas in contrast, the lipid bilayer restricts the embedded domains. Hence, the energies are even sensitive to the phase of the lipid. The folding process of the R-helix in a membrane follows a two-stage mechanism.23-31 First, the hydrophobic polypeptide chain forms a stable transmembrane R-helix crossing the lipid
10.1021/jp810772a CCC: $40.75 2009 American Chemical Society Published on Web 03/19/2009
Hydrogen Bonds in Membrane Proteins bilayer region. Then, these helical domains self-assemble laterally to form the native transmembrane protein. For a typical transmembrane helix, the hydrogen bonds are intrahelical, with interactions between residue i and residue i + 4, although some deformed R-helical structures with hydrogen bonds between residue i and residue i + 5, called the π-bulge, also exist. Some typical examples of a π-bulge are revealed in light-driven anionpumped halorhodopsin and light-driven proton-pumped bacteriorhodopsin.32,33 This distortion leads to a local conformational instability, which creates a hot spot in the protein for conformational change and disruption in the hydrogen-bonding pattern of the R-helix. The lipid bilayer acts to restrict some of the structural diversity of transmembrane domains in membrane proteins. Although, typically, the hydrogen bond dominates the stability and specificity of protein structure, the burial of a polar side chain in an aqueous environment is not always compensated by hydrogen-bond formation; that is, the hydrogen bond cannot always provide enough energy for the protein structural formation. Thus, polar interactions contribute more to specificity than to stability in protein formation. However, on the other hand, in nonpolar environments such as the biological membrane system, hydrogen bonds stabilize the helical conformation of the membrane protein domain and its tertiary interactions.34-36 The typical energetics of NsH · · · O hydrogen bonds is as much as 4.0-5.5 kcal/mol in vacuum, whereas for a solvated protein, because of the extrusion of the polar side chain into the bulk aqueous environment and the solvation of the hydrogen bond by the bulk medium, the hydrogen-bond energy is lowered to some 1.5 kcal/mol.37 In an apolar biological membrane with a dielectric constant close to ε ) 1, the membrane-protein hydrogen-bond energy should be closer to that of the hydrogen bond in vacuum. Theoretical studies of the hydrogen bond in the gas phase have mostly focused on quantum computations for small molecules.38 Various methods have been used such as density functional theory (DFT) and Hartree-Fock methods. For protein molecules, DFT has been applied to study only the side-chains parts,39,40 but those ab initio computation methods started from a static configuration and calculate an energy difference by varying the hydrogen-bond distance and structure angles. According to this quantum mechanical computation method, the hydrogen-bond energy depends not only on the bond distance but also on the bond angle and is better fitted by a Morse function.41 Recently, Mills and Dean surveyed the crystal structures in Cambridge Structure Database and found an interesting angle dependence of the hydrogen bond.42 The truncation or renormalization of the hydrogen-bond energy from the gas phase to the solvated environment is extremely large because the hydrogen-bond distortion in the transition state penetrates deep into water. Experimentally, to investigate protein hydrogen-bond formation, a kinetic isotope method has been used.43 Hence, understanding the hydrogenbond energy other than in the gas phase requires the consideration of a nonstatic method. Because large numbers of molecules are involved in solvated systems, DFT computations are not feasible. This difficulty is aggravated by the fact that one must consider the same calculation for the activated complex. These facts have motivated us to adopt the molecular dynamics simulation approach, although the force field used in the MD simulations is fitted by a gas-phase ab initio computation. The kinetic method contains the ensemble-averaged information. Therefore, the angular part is averaged out in this method.
J. Phys. Chem. B, Vol. 113, No. 15, 2009 5319 The situation for lipids is somewhat different. Each lipid molecule consists of two major parts: One end is the polar head, and the opposite end is a nonpolar tail. The nonpolar tails interact with each other to form a bilayer structure, and the polar heads are pointed toward the solvated water system. Also, the nonpolar tails look like branched chains of two linear chains and form phase structures in the lipid bilayer with three different tilt angles. The dielectric environment in the membrane is similar to that in the vacuum. The biological membrane, however, also has a pseudophase transition close to room temperature because of the tilt change of the lipid chain axis. This constitutes a particular challenge to this method. Essentially, an elongated fatty acid has a length of ca. 20 Å for C18. Once the fatty acid is embedded in a lipid bilayer, it would have a thickness of ca. 40 Å. Furthermore, its polar head would cover a lipid area of ca. 64 Å2. It is well known that, for crystalline alkanes, the chain molecule is completely extended and packed together. When the temperature is increased, a sequence of solid-solid phase transitions below the crystal melting point is observed. In each of these phases, the hydrocarbon chains permit a higher mobility, i.e., rotate more freely. In the lipid bilayer, similar phase transitions occur. At low temperature, all of the molecular chains are highly ordered and referred to as the LB phase.44 Increasing the temperature, a first detectable subtransition temperature, TS, is reached, and the subgel phase LC is obtained. A further increase in temperature leads to a second pretransition temperature, TP. There, one then obtains the gel phase Lβ where the polar head groups are disordered. Continuing to increase the system temperature, one reaches the real melting temperature, Tc, which is a sharp transition. Above Tc, the lipid bilayer is in the lamellar liquid crystalline phase, i.e., LR state. In this molten state, the lipid bilayer structure still stays intact, and the fatty acid chains undergo nearly free twisting motions. At low temperatures, these nonpolar tails constitute a crystal structure. With increasing system temperature, the crystal structure melts to form a gel-like structure. Above the phase transition temperature Tc, the nonpolar tail part turns into a fluid structure with a more random configuration. In our present work, we investigate the hydrogen-bond energetics and even demonstrate their sensitivity to these phase transitions. This again demonstrates the importance of the environments. In a biological system, the membrane participates in many vital transport processes such as the permeability of neutral molecules. Above the melting temperature Tc, fatty acid chains are free to rotate, and kinks are formed that enable small molecules to jump to the voids created. The free motion of lipid molecules also helps the lateral diffusion of lipids in bilayers and proteins on cell surfaces. Therefore, the lipid molecule must be mobile, which is facilitated by the lipid structure. The fatty acid chain structure can shift the melting temperature Tc and, hence, depress its own motion. For example, the hydrogen bonds between the head groups of phosphatidylethanolamine actually stabilize the membrane and elevate Tc. Furthermore, lowering the system pH value has the same effect. In the present work, we use molecular dynamics (MD) simulations to directly observe the dynamics of the protein hydrogen-bond rupture process. Our goal here is to extend the kinetic theory developed in our previous work37 to yield the hydrogen-bond energy for the membrane protein. We even found a first-order phase transition in our membrane systems displayed by our MD simulation in which Tc is ca. 290 K, a value very close to the experimental result. We also compare the hydrogen-bond energies above and below Tc. We found that
5320 J. Phys. Chem. B, Vol. 113, No. 15, 2009 the hydrogen-bond energy above Tc is slightly lower than the hydrogen-bond energy below Tc. Above Tc, lipids behave in a more fluidlike manner, and the hydrogen-bond energy is close to the behavior in the water. However, below Tc, the rigid lipid molecules block the motion of the protein and thus prevent longrange interactions, thereby raising the hydrogen-bond energy. Note also that there is little contribution from the entropy here. In this effort to elucidate the behavior of hydrogen bonds in biological membrane systems, we describe our simulation method and a simple model to extract the total hydrogen-bond free energy in section. We then discuss the various environments and the different phases depending on temperatures in section. Finally, conclusions are given in section. II. Simulation Methods Model of the Membrane System. To build up the lipid-protein system, we followed the method of Roux et al.45-47 Initially, we generated a system consisting of effective lipid spheres with a cross-sectional area of 64 Å2 and placed these spheres on planes within the cell at 7 and 17 Å for the DMPC system or 9 and 19 Å for the DPPC system. Using molecular dynamics and energy minimization, randomized final sphere positions were obtained. Then, the spheres were replaced by DMPC (dimyristoylphosphatidylcholine) or DPPC (dipalmitoylphosphatidylcholine) molecules, which were randomly chosen from a preequilibrated DMPC or DPPC coordinate database, and were preequilibrated and prehydrated by 20 water molecules. The DMPC center of mass was placed at the sphere center. This was followed by a rigid rotation of the lipid chain along the z axis and translation in the xy plane to ensure that bad contacts were removed and that the energy was minimized. Then, the lipid-protein system was hydrated by pre-equilibrated water molecules, and energy minimization was used to refine the configurations. Particular care was required during the placement of each lipid to avoid any bad van der Waals contacts that could give rise to very large gradients/forces during the minimization and equilibration phases. Each monolayer was composed of 48 DMPC lipid molecules with a surface area of ∼64 Å2/lipid. Computational Procedures. To investigate the hydrogenbond rupture process of proteins embedded in the membrane system, molecular dynamics simulations were performed using the program CHARMM48 and the CHARMM27 force field49 with the TIP3P water model.50 Two initial structures of melittin and truncated-KcsA proteins with R-helix structures were taken from the Protein Data Bank with codes 2mlt and 1bl8, respectively. Melittin, with 26 residues, is the main component in bee venom and is a hemolytic and amphiphile peptide. KcsA is a K+ channel membrane protein from Streptomyces liVidans, and in this work, KcsA was truncated from the sequence Val75Ala-Gly-Ile-Thr-Ser-Phe-Gly-Leu- Val-Thr-Ala-Ala-Leu-AlaThr-Trp-Phe-Val-Gly94, named KcsA20. One β-hairpin structure with 12 residues (12-mer), i.e., Ace-(V)5PG(V)5-NH2 with Ace ) sCOCH3, was built. Each peptide was executed by well minimization. The system was built by inserting the minimized peptide into 48 DMPC- or DPPC-based lipid bilayers, including around 1849 or 7373 water molecules, respectively. Because the melittin peptide was horizontally laid on the xy plane and embedded in the DMPC bilayer, without fixing the melittin position, it would move toward the water-lipid interface with the segment Lys21-Arg-Lys-Arg-Gly-Gln26 exposed to water and the rest retained inside the lipid. To illustrate the orientation effect, β-hairpin and KcsA20 were inserted perpendicularly to the membrane surface in the DMPC and DPPC lipid bilayer.
Sheu et al. Periodic boundary conditions were applied to simulate a planar layer in the x and y directions and a bilayer in the z direction with a cell of 32 × 48 × 60 Å3 for the DMPC system or 60 × 60 × 100 Å3 for the DPPC system. The protonation states of ionizable residues were set for pH 7. No counterions were added to the system. The nonbonded interaction was truncated at 12 Å, and the particle mesh Ewald (PME) method51 was employed to treat the long-range electrostatics and to reduce the negative effects of a cutoff. The nonbonding pair list was updated every 50 steps. The SHAKE algorithm52 was applied to constrain the bond lengths involving hydrogen atoms. A time step of 1 fs was used to integrate the Langevin equation. In the lipid system, the simulation cell was anisotropic, and the MD simulations were constrained in the NPT ensemble (i.e., constant pressure and temperature). All of the simulations were performed at 250, 285, 300, and 350 K. The gel-liquid phase transition temperature (Tc) is 296 K for DMPC and 303 K for DPPC,53,54 so two of the temperatures chosen are above Tc and the other two are below Tc in order to ensure that we can study the activation processes in two different phases. Each lipid-peptide system was separately subjected to a 2000-step CONJ energy minimization, followed by 2 ps of molecular dynamics to heat the system temperature to the target temperature. Subsequently, an initial 125-ps MD simulation was performed for equilibration of the system, and the system was coupled to a heat bath with a fixed temperature using Langevin dynamics. A harmonic constraint was used to prevent largescale motion of the protein and to keep the planar shape of the lipid structure. Also, water penetration into the lipid bilayer was inhibited in this way. The harmonic constraint was gradually removed to ensure that a completely free system was achieved within the 50-ps MD simulation. In the following time steps, we allowed the system velocity to be rescaled to adjust the system temperature. A further 5-ns MD equilibration was performed to equilibrate the lipid-peptide system. After the system was equilibrated, we performed MD simulations for an additional 125 or 500 ps and saved the trajectory every 10 fs to ensure that we had enough time resolution for the analysis of the hydrogen-bond rupture process. A time series of the HN · · · O distance of each hydrogen-bond pair was analyzed during the MD simulations. III. Kinetic Model Typically, the standard method for estimating the activation energy is to sum over all thermodynamic structures. The average is taken as a sum over all bond structures, that is, an ensemble average. This method can be referred to as the thermodynamics method, which then must be applied to the protein solvation process, in which there are many bond rupture and association processes. In some situations, it would be better to look at individual hydrogen-bond rupture processes rather than to sum all of them together. Hence, we examine here the detailed microscopic processes in protein dynamics and count each individual process in a molecular dynamics experiment as a direct input from rate theory as introduced in our previous work.37 A classical Kramer’s theory estimates the frictiondependent transition rate by counting a successful trajectory as climbing from the well minimum to reaching the barrier top. In this way, we establish a connection between the reaction rates and mean first passage time (mft), and we can then extract the activation energy through a computation of the mean first passage. Based on this method, which we term a kinetic method, one can observe, through a time series, the microscopic hydrogen-bond rupture and formation processes.
Hydrogen Bonds in Membrane Proteins
J. Phys. Chem. B, Vol. 113, No. 15, 2009 5321
Figure 1. Schematic potential surface for the protein hydrogen-bond rupture process.
Protein hydrogen-bond rupture processes can be imagined as a unitary reaction with the dissociation energy curve shown in Figure 1. At the minimum point A, the hydrogen bond stays at its equilibrium state, and its typical distance is ca. 2.0 Å. Thermal motion keeps the peptide chain fluctuating until the hydrogen-bond distance is increased to 3.5 Å, i.e., the barrier top B. This entire process is counted as a rupture process for a solvated system. It is similar to a Brownian particle climbing a barrier top. Our argument is based on Yamamoto’s correlation function expression, i.e., a microscopic average of the fluctuation-dissipation theorem
k(t) )
〈hA[x(0)] · hB[x(t)]〉 ≈ kAfB exp(-t/trxn) 〈hA[x(0)]〉
(1)
Here, hA[x(0)] is a characteristic function that counts the number of trajectories reaching x(0) at point A, and hB[x(t)] in the same way counts the number of trajectories reaching x(t) at point B. Using a random trajectory with a time scale longer than the system relaxation time trxn, one can separate the factors of the rate constant kAfB and the relaxation process from the timedependent rate constant k(t). In this convention, x(t) represents a point in phase space at a certain instant, and 〈 · · · 〉 denotes the ensemble average in phase space. Two remarks are worth mentioning. First, here, our molecular dynamics simulation is based on Yamamoto’s picture, so the measured rate constant is a microscopic ensemble average of all of the degrees of freedom of hydrogen-bond rupture, that is, an ensemble average of all types of H-bond motions. Our estimation of the HN · · · O distance actually averages out the bending of the NH · · · O angle. Next, our MD simulations are performed at a constant pressure and constant temperature. In our simulations, the hydrogen-bond length cutoff is 3.5 Å for the N · · · O distance. IV. Activation Thermodynamics Based on the preceding Yamamoto scheme, one is able to count the rupture time and then to estimate the rate. In this argument, we are then ready to fit the activation thermodynamic quantities. Here, we briefly review the connection. There is also in the solution an entropy effect because of the large-scale motion during the rupture processes. Chemical reaction in solution depends on many phenomenological parameters such as temperature, friction and pressure etc. Experimentally, one can then extract the thermodynamic + energy from the Arrhenius equationk ) Ae-Ea /kBT. Here, the prefactor A is essentially temperature-independent. The activa-
tion energyE+ a , which is the internal energy needed for breaking the bond in the reaction, dominates the temperature-dependent effect. kB is the Boltzmann constant. According to the measured reaction rates, A and E+ a are obtained by plotting log k versus 1/T (Arrhenius plot), from which various kinds of activation thermodynamics quantities can be obtained. We first explore a gas-phase rate theory in Figure 1 as an example where E+a indicates the barrier height for escaping from well A to the dissociation state, where the basic assumption here is that the system is in equilibrium at a temperature T during the entire reaction. The molecule moves in well A, where it is initially located, with a frequency given by kBT/h, which is about 1013 s-1 at 300 K. According to transition state theory, the chemical reaction rate is a product of an oscillation frequency, the degeneracy of the state at the barrier top, and the activation + process, i.e., kTST ) νZe-Ea /kBT, where Z is the ratio of the number of states or the degeneracy of the dissociation state to the number + of states in well A and is equal to eSa /kB, where S+ a is the entropy difference between the well and the dissociation state. We then + + + have kTST ) νe-Fa /kBT, where F+ a ) Ea - TSa is the difference in the Helmholtz energy between well A and the dissociation state and is called the activation Helmholtz energy. The experiment is done under a pressure P and the volume of the system is denoted as VA in state A and VAB in the dissociation state. The particle climbing from well A to the dissociation state + + requires an energy of E+ a + PVAB, with the volume change VAB + ) VAB - VB, in order to reach the dissociation state, where VAB is called the activation volume. We can summarize this analysis + + ) E+ to obtain the activation enthalpy HAB a + PVAB. The simplest expression of the gas-phase transition-state rate constant + is then kTST ) νe-GAB/kBT, with the activation Gibbs free energy + + + . According to the Gibbs free energy GAB ) HAB - TSAB argument, we reassemble the activation entropy to the prefactor + part A ) νe-SAB/kB. Using ν ) kBT/h, we obtain the Eyring + expression kTST ) (kBT/h)e-GAB/kBT, where the activation part is determined by the Gibbs free energy rather than by the internal energy as in the Arrhenius expression. For the simplest example in a gas-phase reaction, PV ) nRT, where n is the number of moles, R is the gas constant, and the PV difference is defined between the well A and the dissociation state as ∆(PV) ) ∆nRT, where ∆n is the difference in the number of moles for the reaction process. For the dissociation reaction C f A + B, then ∆n is equal to 1. We finally have the corresponding activation + + energies with the Arrhenius law E+ a ) HAB - RT and SAB ) R ln A - R ln T - R ln(kB/h). Now, we turn to consider a chemical reaction in solution, which has a more complicated behavior. A phenomenological theory is used to adopt macroscopic parameters such as temperature, friction, and so on. This simplified formulation based on statistical mechanics can provide an intrinsic understanding of the chemical reaction. The easiest way to include the solvent dynamic effect is to adopt Kramer’s theory. A phenomenological friction is included in this theory. To introduce Kramer’s theory, we first introduce some conventions. Here, the nonequilibrium statistical mechanics relationship between the velocity autocorrelation time τv and the friction coefficient ζ is established as τv ) m/ζ, where m is the particle mass. Based on the Stokes law in fluid mechanics, there is a simple relationship between viscosity η and friction, namely, ζ ) 6πηa, where a is the particle radius. The effect of friction is illuminated and characterized by the important transmission coefficient κ and k ) κkTST, where κ includes the solvent dynamic effects other than the gas-phase transition-state rate information. Comparing the rate constant with friction to the
5322 J. Phys. Chem. B, Vol. 113, No. 15, 2009
Sheu et al.
TABLE 1: Mean First Passage Times (mft’s) and Activation Energies of the Peptides at Various Temperatures temperature (K)
mft (ps)
+ HAB (kcal/mol)
prefactor (s-1)
Melittin,a R-Helix in DMPC,b H-Bond Pair Ile2-O · · · HN-Leu6 (Horizontal Orientation) 250 59.43 3.61 2.39 × 1013 285 24.37 300 19.23 2.62 4.20 × 1012 350 10.27 KcsA20,c R-Helix in DPPC,b H-Bond Pair Leu27-O · · · HN-Arg30 (Vertical Orientation) 250 78.85 3.80 2.68 × 1013 285 30.80 300 23.78 2.67 3.70 × 1012 350 12.54 β-Hairpind in DMPC, H-Bond Pairs Val10-O · · · HN-Val4 and Val11-O · · · HN-Val5 (Vertical Orientation) 250 62.0 3.48 1.77 × 1013 285 26.23 300 21.0 2.68 4.26 × 1012 350 11.05 a Melittin: PDB code 2mlt. b DMPC, dimyristoylphosphatidylcholine; DPPC, dipalmitoylphosphatidylcholine. The gel-liquid phase transition temperature is 296 K for DMPC and 303 K for DPPC.53,54 c KcsA20: structure truncated from PDB code 1bl8 with the sequence Val75-Ala-Gly-Ile-Thr-Ser-Phe- Gly-Leud Val-Thr-Ala-Ala-Leu-Ala-Thr-Trp-Phe-Val-Gly94. β-hairpin: 12-mer residues with the sequence Ace-(V)5PG(V)5-NH2.
+
Arrhenius equation, we obtain the connection k ) κνe-GAB/kBT + ) Ae-EAB/kBT. Because E+ a is the internal energy, the prefactor for a chemical reaction in solution is then equal to A ) κν + + exp[SAB /kB - ∆(PVAB )/kBT]. For a unitary reaction, A f B, the activation volume is almost fixed; hence, we can ap+ ) ≈ 0. We therefore obtain the standard proximate ∆(PVAB + + thermodynamic quantities as E+ a ≈ HAB and SAB ) R ln A - R ln κ - R ln T - R ln(kB/h). This suggests a convenient way to obtain the activation energy and the prefactor through an + Arrhenius fit; the accompanying quantities of interest are HAB + and SAB. The Gibbs free energy then satisfies the relationship + + + ) HAB - TSAB . GAB Accordingly, the entropy effect for the unitary reaction is included in the prefactor part. A typical prefactor value is ca. 1013 s-1; hence, the solvent dynamics effect and a decreasing of the entropy effect can suppress the A value. For the unitary reaction, our H-bond energy is equal to the enthalpy, as well. V. Results and Discussion The most commonly used lipid molecules in molecular dynamics simulations are DPPC and DMPC, among others. There exist different phases for lipids that have fluidlike and gel-like structures. This gel-fluid phase transition temperature, Tc, is close to room temperature for the lipid molecules of interest (see Table 1). The hydrogen-bond energy is defined as the free energy or the average of energy fluctuations. Different phases of the lipid membrane exhibit different protein hydrogenbond length fluctuations, although this is not a true critical phenomenon and is sometimes referred to as a pseudophase transition. Our molecular dynamics simulation was done in different temperature regimes, namely, above Tc and below Tc. We expect that, in the fluid phase, the lipid molecule should have more freedom to move (i.e., twist motion) and there should be larger-scale fluctuations of the rupture process. On the contrary, in the gel-like phase, the lipid molecules cannot move
Figure 2. Structure of melittin in the membrane. The melittin was inserted into the center of the membrane and was oriented parallel to the membrane surface. Melittin is a very short transmembrane protein and has an R-helix structure. We placed this R-helix at the interface between the lipid bilayers of the DMPC membrane.
as freely, and the lipid molecules can block the structure change of the membrane protein. In our previous work,37 we developed a kinetic method to investigate an individual pair of hydrogen bonds through its own rupture process. We found that the hydrogen-bond energy is greatly reduced from ca. 5 kcal/mol in the gas phase to ca. 1.5 kcal/mol in water. Because of the entropy fluctuation of the protein chain ends, the R-helix and β-sheet exhibit different hydrogen-bond energies. This type of entropy fluctuation also exists in the membrane protein system. In our simulations, we used the CHARMM force field TIP3P, which was fitted by ab initio computation. A typical ab initio computation result of the hydrogen-bond energy for the water dimer is ca. 5.0 kcal/mol.55 Melittin is a very short transmembrane protein and has an R-helix structure. We put this R-helix at the interface between the lipid bilayers of the DMPC membrane (see Figure 2). To study the β-sheet structure, we used a β-hairpin as in ref 37 and inserted it into the DMPC membrane (see Figure 3). A typical rupture process is shown in Figure 4, and the hydrogen pairs chosen are listed in Table 2. The two structures have similar hydrogen-bond energies. Before discussing the H-bond energies above and below the membrane transition temperature, we show that a first-order phase transition is directly observed in the calculations of our membrane system. Here again, the same bond is being broken but in two different lipid environments. In Figure 5, we plot the Gibbs free energy vs T, where the free energy is taken from the logarithm of the product of the inverse of the mft and the temperature in Table 3 for KcsA20. We observe an abrupt change of the free energy at ca. 290 K. Note that KcsA20 is observed as an R-helix embedded in DPPC with a Tc of ca. 303 K. Note
Hydrogen Bonds in Membrane Proteins
J. Phys. Chem. B, Vol. 113, No. 15, 2009 5323 TABLE 2: Hydrogen-Bond Pairs for Melittin and β-Hairpin type of interactiona
H-bond pair R-Helix Ile2-O · · · HN-Leu6 Val10-O · · · HN-Val4 Leu21-O · · · HN-Gln26 Gly1-O · · · HN-Val5
O-O O-O O-I I-O β-Sheet
Val10-O · · · HN-Val4 Leu16-O · · · HN-Ile20 Val5-O · · · HN-Leu9 a
Figure 3. Structure of the β-hairpin in the membrane. We compared the effect of the orientation for the R-helix in the membrane and found that there was no difference. Hence, we chose only one orientation for the β-sheet in the membrane.
O-O O-O O-O
O, hydrophobic; I, hydrophilic.
inhomogeneous system, there might be a minor difference in the H-bond energy from pair to pair or even from site to site or + is higher than that in from protein to protein. The enthalpy HAB a water system, which is typically ca. 1.5 kcal/mol. This strong reduction of the strength of the hydrogen bond in water is attributed to the dynamic similarity with the water medium. This reduction results from the interaction of the protein hydrogen bond with the similar long-range energetics in water. This similarity will be reduced in the lipid environment, hence leading to a rise in the activation energy. Presumably, the membrane environment has a dielectric constant close to that of a vacuum, and the rupture process there should have a higher energy, as reported in ref 37, closer to the large hydrogen-bond energy in a vacuum. The entropy term is similar to the vacuum value for the rigid structure below Tc. The reaction rate is low because of the rigidity (Table 3). Below Tc, the rate is high, even higher than in water. This is due to the absence of the entropy bottleneck. Hence, the rates are higher despite a small increase in activation energy. Hence, we have the highest rate constant in lipids, despite the higher activation energy. Hence, the hydrogen-bond chemistry is much enhanced in lipids above Tc, even when compared to the extremely small activation energy in water. This makes lipids into such an attractive medium for biochemical reactions, even when compared to water. Hence, if we compare the hydrogen-bond energies of R-helix and β-hairpin in diverse environments (see Table 4 and the
Figure 4. Time series of the rupture processes for the H-bond pair Leu27-O · · · HN-Arg30 of KcsA20.
that our simulated Tc is very close to the experimental Tc value within an error of ∼4%. This represents an interesting confirmation of this dynamic method of determining energies that are so strongly dependent on the solvent; in fact, it even includes energy changes due to phase transitions. These features are usually not included in the static method of calculating statistical thermodynamic quantities of reactants and products. In our simulations above Tc, both the R-helix and the β-sheet (see Table 1) show larger hydrogen-bond energies (enthalpies) than the results in water (see Table 4). Note that, in Table 4, we report the hydrophobicities of the H-bond pairs. Most of them are hydrophobic pairs. Because a protein is a highly
Figure 5. Plot of the Gibbs free energy (in units of the Boltzmann constant kB) vs temperature. The Gibbs free energy is obtained as the logarithm of the product of the inverse of the mean first passage time (mft) and the temperature. The mft is taken from Table 3 for KcsA20. An abrupt change of the free energy occurs at ca. 290 K, which shows a phase transition temperature, Tc. Note that the KcsA20 is an R-helix embedded in DPPC with Tc ≈ 303 K. Our simulated Tc value is very close to the experimental Tc value with an error of ∼4%.
5324 J. Phys. Chem. B, Vol. 113, No. 15, 2009
Sheu et al.
TABLE 3: Hydrogen-Bond Rupture Times (ps) in vacuum37 R-helix
temperature (K)
R-helix
β-hairpin
250
300
221.32
350
63.9
104
400
22.32
44
in H2O37
R-helix
β-hairpin
R-helix
21
85.3
100.3
46
40
11.046
35.9
48
29
27.4
β-hairpin
78.85 (59.43) 30.80 (24.37) 23.78 (19.23) 12.54 (10.27)
285
a
in D2O22
in lipid a
β-hairpin
62 26.23
Data in parentheses are just from a different run.
TABLE 4: Prefactors (s-1) and Activation Energies (kcal/mol) of the Hydrogen-Bond Rupture Processes in Various Phases in vacuum37
in D2O22
in lipid
+ HBA
A
+ HBA
A
5.55
5.02 × 1013
3.80a,b 3.61b,c 2.67d 2.62d
2.68 × 1013 2.39 × 1013 3.70 × 1012 4.20 × 1012
4.79
9.37 × 1012
3.48b 2.68d
in H2O37
+ HBA
A
+ HBA
A
3.61
5.01 × 1012
1.93
5.49 × 1011
β-Hairpin 1.77 × 1013 3.07 4.26 × 1012
1.73 × 1012
1.58
3.53 × 1011
R-Helix
a
Melittin. b Below Tc. c KcsA20. d Above Tc.
associated rupture times in Table 3). Our finding of the hydrogen-bond energy in a lipid is lower than that in a vacuum but higher than that in water. The hydrogen-bond energies follow the order vacuum > lipid > D2O > H2O. However, the reduction of the A factor in water is nearly 2 orders of magnitude, whereas in lipids it is near unity. Orientation Effect. It would be quite interesting to know whether the hydrogen-bond strength is orientation-dependent, that is, dependent on whether the transmembrane protein is perpendicular or parallel to the lipid molecule. In our R-helix simulation, we allowed the melittin to be parallel to the membrane. We then adopted another transmembrane protein, KcsA20, with one of its bundles perpendicular to the membrane. The results are presented in Table 1. Comparing the results for the R-helix and KcsA20, one van see that both the prefactor and + values are quite similar. Therefore, the orientation of the HAB protein bundle does not influence the H-bond. VI. Conclusions The protein hydrogen bond is a key element in defining protein structure and perhaps more important in defining different dynamics of protein motion in various environments, including the motions of folding. One of the important properties of such hydrogen bonds seen here appears to be that their properties are only definable holistically as a unit together with their environment. This environment changes the chemistry of these hydrogen bonds beyond the effects predicted by molecular partition functions. As an example, in ref 37, we showed the fundamental difference between protein hydrogen bonds in the isolated state and in water systems. There is a substantial reduction of the H-bond energy from its inherent value in a vacuum of some 5 kcal/mol to a much smaller value in a water environment of some 1.6 kcal/mol. Thus, there is a very strong solvent effect on the strength of a hydrogen bond. This finding has been supported recently by molecular dynamics measurements in D2O yielding a very large solvent isotope effect. More
fundamental and novel is the direct observation of a very large entropy effect in water as a result of the large dynamically clustered water environment. This even leads to an entropy bottleneck reducing the rate in water by some 2 orders of magnitude despite the very low activation energy. There is no such bottleneck for the isolated molecule. This is the result of a very large negative change in the entropy of activation in water. Therefore, the hydrogen bond must be taken as a unit together with this particular environment. This entropy cannot be directly obtained from the partition functions of the molecules involved, but rather is due to large long-range environmental effectssvery often seen here to be of dominant magnitude for the ensuing chemistry. There is a large increase in rate due to the drop in energy but at the same time a 100-fold entropic bottleneck opposing it. This entropy carries with it interesting possibilities of modification with further external agents or pH in water. The low energy in this case is very close to that postulated for the energetics deemed operative in the fluctuation equilibrium in pure water, a process that might even be ratedetermining in this case and explains the ready energy transfer from the solvent medium. Switching to the lipid environment creates a low-energy path for carrying out protein chemistry, such as in complex enzymes, of some 3 kcal/mol of binding energy as required for enzymes. In the heterogeneous membrane environment, the intraprotein hydrogen bond can turn its direction toward the lipid molecule, even though the dielectric constant is only slightly different from that of the vacuum system. Lipid molecules with alkane chains also change the hydrogen-bond direction. The solvation effect also depresses the energy of the hydrogen bond in lipids, although not to the level observed for water. The work carried out here shows that the hydrogen bond in a membrane environment constitutes yet a third class of the cooperation of the hydrogen bond with the medium. We observe that the hydrogen-bond energy is strongly reduced to some 3.6 kcal/mol in the membrane medium. In the fluid phase above
Hydrogen Bonds in Membrane Proteins the transition temperature of the lipid, both the R-helix and the β-sheet have slightly smaller hydrogen-bond energies. Interestingly, the entropy bottleneck now nearly disappears, leading to a near-normal A factor and again a smaller prefactor in the fluidlike phase than in gel-like phase, which has slower motion and less energy fluctuation. The result now seen for lipids is that the activation entropy is closer to that found for the vacuum. This leads to a most interesting consequence. Even though the activation energy in the lipid is higher than that in water, the lack of the entropy bottleneck now enhances the rate above that found for water (Table 3). This rate now is higher than in pure water. Because this energy, however, is comparable to some enzyme complex energies, these chemical processes now become reasonable processes to occur. Hence, substantial chemistry can occur because of the stronger hydrogen bonds, without giving up the reactivity of the hydrogen bond, the latter being required in addition for dynamic protein processes in the chain. This puts lipids into a very special class as a reactive media. This brings us to three prototypes of the behavior of hydrogen bonds acting in concert with their medium to substantially define the chemical properties of the hydrogen bond/complex. The environment of the hydrogen bond must therefore be taken as a unit, and only then can the completely different possible chemistries for these three cases be explained. Therefore, in summary, the three cases are as follows: 1. The isolated molecule with no entropy bottleneck has an energy of some 5 kcal/mol. This is well-documented from ab initio theory and in these MD calculations. This value is too high for reaction to occur at normal physiological temperatures, and hence, the system will remain intact. 2a. In the case of the protein in H2O, a very large long-range fluctuating water environment directly interacts with the hydrogen bond to lower its energy to the well-known low value of 1.6 kcal/mol and thus facilitates ready hydrogen-bond breaking and making within the protein. 2b. The entropy bottleneck for the water complex leads to a reduction in rate of some 2 orders of magnitude. This low value of the energy will, however, still lead to very facile structural rearrangements in the protein. 3. For the protein in a lipid environment, the energy is reduced to some 3.6 kcal/mol, but now there is no entropy bottleneck. Hence, although the energy barrier is higher than in case 2, the lack of a bottleneck makes the rate higher than in case 2. The higher barrier is particularly attractive to enzyme reactions, where the enzyme complex typically has an energy above some 3 kcal/mol, and hence, the kinetic competition between these processes enhances enzyme reactivity. This might be similar to the case of organic solvents, which tend to increase enzyme rates. The weakness of the hydrogen bond demands a holistic picture including the global environment as a unit. Without this environment, hydrogen bonds are essentially rigid. Such a model can even display a correct phase transition. Herein, we discussed prototypes for very large environmental influences on the hydrogen bond in proteins. Any actual case might need further particulars in the model, such as pH, denaturants, and so on. Acknowledgment. S.-Y.S. and D.-Y.Y. gratefully thank the National Science Council of Taiwan for Grants NSC-97-2113M-010-002-MY2 and NSC-97-2113-M-001-021-MY2, respectively. The authors are grateful to the National Center for HighPerformance Computing of Taiwan. This work was supported by the Taiwan/Germany program at the NSC/Deutscher Aka-
J. Phys. Chem. B, Vol. 113, No. 15, 2009 5325 demischer Austauschdienst. A portion of the research described in this article was performed in the Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the Department of Energy Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory. This research was supported by the DFG cluster of excellence Munich Centre for Advanced Photonics. References and Notes (1) Arunan, E. Curr. Sci. 2007, 92, 17–19. (2) Pimentel, G. C.; McClellan, A. L. The Hydrogen Bond; W. H. Freeman & Co.: San Francisco, CA, 1960. (3) Theoretical Treatments of Hydrogen Bonding; Hadzˇi, D., Ed.; John Wiley & Sons, Chichester, U.K., 1997. (4) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press: Oxford, U.K., 1997. (5) Coulson, C. A. Research 1957, 10, 149–159. (6) Buckingham, A. D.; Del Bene, J. E.; McDowell, S. A. C. Chem. Phys. Lett. 2008, 463, 1–10. (7) Mattos, C. Trends Biochem. Sci. 2002, 27, 203–208. (8) Doruker, P; Bahar, I Biophys. J. 1997, 72, 2445–2456. (9) Chakraborty, S.; Bandyopadhyay, S. J. Phys. Chem. B 2007, 111, 7626–7630. (10) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, J. J. Phys. Chem. B 2006, 110, 4393–4398. (11) Patel, S.; Mackerell, A. D, Jr.; Brooks, C. L., III J. Comput. Chem. 2004, 25, 1504–1514. (12) Sheu, S. Y.; Yang, D. Y.; Selzle, H. L.; Schlag, E. W. J. Phys. Chem. A 2002, 106, 9390–9396. (13) Perrin, C. L.; Nielson, J. B. Annu. ReV. Phys. Chem. 1997, 48, 511–544. (14) Rozas, I. Phys. Chem. Chem. Phys. 2007, 9, 2782–2790. (15) Dannenberg, J. J. J. Mol. Struct. (THEOCHEM) 1997, 401, 279– 286. (16) Jeffrey, G. A.; Saenger, W. Hydrogen Bonding in Biological Systems; Springer Verlag: Berlin, 1991. (17) Sheu, S. Y.; Selzle, H. L.; Schlag, E. W.; Yang, D. Y. Chem. Phys. Lett. 2008, 462, 1–5. (18) Klotz, I. M. Protein Sci. 1993, 2, 1992–1999. (19) Williams, D. H.; Searle, M. S.; Mackay, J. P.; Gerhard, U.; Maplestone, R. A. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 1172–1178. (20) Fersht, A. R.; Shi, J.-P.; Knill-Jones, J.; Lowe, D. M.; Wilkinson, A. J.; Blow, D. M.; Brick, P.; Carter, P.; Waye, M. M. Y.; Winter, G. Nature 1985, 314, 235–238. (21) Davis, A. M.; Teague, S. J. Angew. Chem., Int. Ed. 1999, 38, 736– 749. (22) Sheu, S. Y.; Schlag, E. W.; Selzle, H. L.; Yang, D. Y. J. Phys. Chem. A 2008, 112, 797–802. (23) Ubarretxena-Belandia, I.; Engelman, D. M. Curr. Opin. Struct. Biol. 2001, 11, 370–3076. (24) Senes, A.; Ubarretxena-Belandia, I.; Engelman, D. M. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 9056–9061. (25) Killian, J. A.; Nyholm, T. M. Curr. Opin. Struct. Biol. 2006, 16, 473–479. (26) Ja¨hnig, F. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 3691–3695. (27) Toraya, S; Nagao, T; Norisada, K; Tuzi, S; Saito, H; Izumi, S; Naito, A. Biophys. J. 2005, 89, 3214–3222. (28) Tierney, K. J; Block, D. E; Longo, M. L. Biophys. J. 2005, 89, 2481–2493. (29) White, S. H.; Wimley, W. C. Annu. ReV. Biophys. Biomol. Struct. 1999, 28, 319–365. (30) Kim, S; Cross, T. A. Biophys. J. 2002, 83, 2084–2095. (31) Chakrabarti, P.; Chakrabarti, S. J. Mol. Biol. 1998, 284, 867–873. (32) Lanyi, J. K. J. Struct. Biol. 1998, 124, 164–178. (33) Oesterhelt, D. Angew. Chem., Int. Ed. Engl. 1976, 15, 17–24. (34) White, S. H. AdV. Protein Chem. 2005, 72, 157–172. (35) Grigoryan, G.; DeGrado, W. F. Nat. Chem. Biol. 2008, 4, 393– 394. (36) Joh, N. H.; Min, A.; Faham, S.; Whitelegge, J. P.; Yang, D.; Woods, V. L, Jr.; Bowie, J. U. Nature 2008, 453, 1266–1270. (37) Sheu, S. Y.; Yang, D. Y.; Selzle, H. L.; Schlag, E. W. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 12683–12687. (38) Kim, K.; Friesner, R. A. J. Am. Chem. Soc. 1997, 119, 12952– 12961. (39) Morozov, A. V.; Kortemme, T.; Tsemekhman, K.; Baker, D. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 6946–6951. (40) Li, J. L.; Car, R.; Tang, C.; Wingreen, N. S. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 2626–2630. (41) Kang, Y. K. J. Phys. Chem. B 2000, 104, 8321–8326.
5326 J. Phys. Chem. B, Vol. 113, No. 15, 2009 (42) Mills, J. E. J.; Dean, P. M. J. Comput.-Aided Mol. Des. 1996, 10, 607–622. (43) Krantz, B. A.; Srivastava, A. K.; Nauli, S.; Baker, D.; Sauer, R. T.; Sosnick, T. R. Nat. Struct. Biol. 2002, 9, 458–463. (44) Metzler, D. E. Biochemistry: the Chemical Reactions of LiVing Cells; Academic Press: San Diego, CA, 2001. (45) Woolf, T. B.; Roux, B. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 11631–11635. (46) Woolf, T. B.; Roux, B. Protein Struct. Funct. Genet. 1996, 24, 92–114. (47) Berneche, S.; Nina, M.; Roux, B. Biophys. J. 1998, 75, 1603–1618. (48) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminarhan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187–217. (49) MacKerell, A. D., Jr; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.;
Sheu et al. Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586–3616. (50) Jorgensen, L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1993, 79, 926–935. (51) Darden, T.; York, D.; Pederson, L. J. Chem. Phys. 1993, 98, 10089– 10092. (52) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (53) Lee, C.-H.; Lin, W.-C.; Wang, J. Phys. ReV. E 2001, 64, 20901– 20904. (54) Taly, A.; Baciou, L.; Sebban, P. FEBS Lett. 2002, 532, 91–96. (55) Scheiner, S. Annu. ReV. Phys. Chem. 1994, 45, 23–56.
JP810772A
