Free-energy profiles describing the relative orientation of membrane proteins along predefined coordinates can be efficiently calculated by means of umbrella simulations. Such simulations generate reliable orientational distributions but are difficult to converge because of the very long equilibration times of the solvent and the lipid bilayer in explicit representation. Two implicit lipid membrane models are here applied in combination with the umbrella sampling strategy to the simulation of the transmembrane (TM) helical segment from virus protein U (Vpu). The models are used to study both orientation and energetics of this α-helical peptide as a function of hydrophobic mismatch. We observe that increasing the degree of positive hydrophobic mismatch increased the tilt angle of Vpu. These findings agree well with experimental data and as such validate the solvation models used in this study.
When carrying out molecular dynamics simulations of biological systems with classical force fields (Karplus and McCammon 2002), the most accurate representation involves modeling both the macromolecules and the solvent explicitly at atomic-level detail. However, because the solvent molecules dominate in number, the majority of the calculation time is spent on solvent–solvent interactions. This therefore limits the feasible amount of conformational samples of the macromolecular solute one is interested in. The problem is particularly exacerbated when the macromolecule undergoes large conformational changes that require significant solvent equilibration time. To this end, much effort has been expended to develop implicit solvent models that would represent the important physical information on the solute–solvent interaction as some effective energy function obtained by formally integrating, under various approximating assumptions, over all the solute degrees of freedom (Roux and Simonson 1999; Feig and Brooks 2004; Brannigan et al. 2006). The fact that the solute degrees of freedom have been integrated out enables for an adiabatic (i.e., instantaneous) response of the solvent to any new conformation of the solute. This instantaneous equilibration is especially important when the macromolecule is solvated a lipid membrane environment, whose slow lateral and transversal reorganization requires relaxation on timescales that are very long and would require excessive computing in an explicit lipid representation (Pastor 1994; Forrest and Sansom 2000).
We focus on two implicit solvent models for lipid bilayer systems. The first model is implicit membrane model 1 (IMM1) (Lazaridis 2003), an extension for proteins in lipid bilayers surrounded by an aqueous environment of a previous effective energy function 1 (Lazaridis and Karplus 1999) developed to model proteins in aqueous solvents. IMM1 has been found to give good results with minimal CPU time and as such is quite promising. Specifically, IMM1 has been applied to the study of structure and energetics of TM proteins (Mottamal et al. 2006), insertion of helices into membranes (Lazaridis 2003), and helix orientation in lipid bilayers (Lazaridis 2003). The second model is the generalized Born with a simple switching model (GBSW) (Im et al. 2003). It builds on previous versions of generalized Born models as approximate solutions to the Poisson–Boltzmann equation. Thus far, GBSW has been applied to the study of helix–helix interaction (Im et al. 2003), helix orientation in membranes, the role of the membrane in dictating conformational dynamics of TM proteins (Im et al. 2003), and interfacial folding of membrane proteins (Im and Brooks 2005). When compared to available experiment data, both IMM1 and GBSW agree reasonably well (see references for details).
Another challenging aspect in simulating biological systems is the determination of free energies changes for solvated systems (Levy and Gallicchio 1998). The majority of methods developed for this purpose are based on thermodynamic perturbation (Zwanzig 1954), thermodynamic integration (Kirkwood 1935), probability distribution methods (Bennett 1976; Torrie and Valleau 1977), or, more recently, the Jarzynski identity (Jarzynski 1997). A particularly useful approach in the category of the probability distribution methods is umbrella sampling (Bartels and Karplus 1997), which allows one to determine the free energies as a function of a predefined reaction coordinate. This is achieved by adding a potential to the Hamiltonian of the system that allows for the uniform sampling of conformational space as parametrized by the reaction coordinate. This methods has been successfully used in many cases, but the process can be time consuming. In cases when the reaction coordinate involves significant conformational displacement in a lipid environment, the ability to generate a uniform sampling during umbrella sampling along that coordinate relies heavily on the ability of the solvent to have sufficient time to respond to the extensive protein change. Such a solvent reorganization is computationally quite demanding, and of particular importance in instances involving the orientation of protein helices inside lipid bilayers. This problem of course does not exist in the case of implicit solvent/lipid models, for which—as a result of the integration over their degrees of freedom—the “solvent” and the “lipid” relax instantaneously in response to any change in orientation or structure.
Hydrophobic mismatch, defined as the difference between the hydrophobic length of TM segments of a protein and the hydrophobic width of the surrounding lipid bilayer, is the driving force behind TM segment reorientation in situations when the hydrophobic stretch of a membrane pass protein is not equal to the thickness of the hydrophobic core of the lipid bilayer (see Killian 1998 for a detailed review). Hydrophobic mismatch can be of two types: positive mismatch, when the hydrophobic stretch of the peptide is longer than the thickness of the hydrophobic core of the membrane, or negative, the case when the thickness of the hydrophobic core exceeds the peptide length. Regardless of its type, this mismatch yields an energetically unfavorable condition: in the case of positive mismatch, some of hydrophobic side chains can get exposed to solvent, and in the case of negative mismatch, some of the hydrophilic sides can find themselves buried within the hydrophobic core of the membrane. Stand-alone TM peptides or TM segments in membrane protein complexes rearrange such as to minimize the hydrophobic mismatch and to achieve hydrophobic matching (White and von Heijne 2005), although in particular instances, interfacial anchor properties of tryptophan residues in transmembrane (TM) peptides can dominate over hydrophobic matching effects in peptide–lipid interactions (de Planque et al. 2003).
In instances of positive mismatch, which is the case of this study, the TM helices rearrange by tilting relative to the membrane normal. Positive hydrophobic mismatching was of interest to us because of the significant amount of experiment data generated about this process. As such, it provided the opportunity to test the applicability of implicit solvent modeling by making a comparison to available experimental data.
The TM helix is the fundamental structural unit of helix bundle membrane proteins. A proper computational description of the free energy of association at the atomic detail is an important ingredient in understanding the details of the amino acid code used by the translocon to select polypeptide segments for insertion as TM helices in helix bundle membrane proteins (White and von Heijne 2005). A good model system for the purpose of our computational study is Vpu, an accessory protein of HIV-1 with well-characterized biological activities; as such, it serves as a suitable subject for the development of methods for structure determination of membrane proteins (Park et al. 2003). Park and Opella (2005) reported liquid and solid-state nuclear magnetic resonance data to derive structural information on the orientation of Vpu. They found that the tilt of a helix in a membrane (in their case the TM protein Vpu) is dependent on positive hydrophobic mismatch. Specifically, they reported that cosine of the tilt angle, θ, increases linearly with the degree of positive hydrophobic mismatch.
In this project, we tested the ability of IMM1 and GBSW to recreate the experimentally observing tilting of the TM helix Vpu in response to positive hydrophobic mismatch. Additionally, umbrella sampling (Bartels and Karplus 1997) was used to obtain free-energy profiles of the tilting as a function of membrane thickness. Assuming that the tilt angle a helix adopts in a membrane represents its energetically favorable orientation, by using these free-energy profiles, we can predict the tilt angle. We can then compare these to the experimental values, and as such evaluate the performance of IMM1 and GBSW solvent models in concert with umbrella sampling. These are critical tests because the tilting requires that the protein systems undergo sweeping motions involving large excursions inside the lipid.
We present the results of this simple study, which show that IMM1 and GBSW are suitable membrane models for carrying out molecular dynamics simulation of simple protein–lipid systems by the fact that umbrella sampling–based free-energy profiles generate hydrophobic mismatch patterns in accord with experimental data.
Methods
The membrane is modeled as a lipid slab between two water regions using two implicit solvent models. The first one, IMM1, uses an empirical energy function, parameterized with experimental data. IMM1 primarily accounts for solvent exclusion effects and assumes that total solvation energy is the sum of group contributions. It partially accounts for electrostatic interaction by using a distance dependent dielectric function allowing for screening of electrical charge. The second model, GBSW, decomposes the solvation free energy into nonpolar and electrostatic components. The nonpolar component is approximated as being proportional to the accessible surface area of the protein. The electrostatic component is determined by an analytical approximation of the Poisson–Boltzmann method by a generalized Born equation for the calculation of the solvation energy for an isolated charge in some dielectric medium.
In our simulations, initial coordinates were obtained from the X-ray structure of TM segment of Vpu (PDB ID: 1pje). For both IMM1 and GBSW, four independent simulations were carried out with a membrane thickness D equal to 15, 19, 23, and 27 Å, respectively. For each simulation, the following protocol was used. The initial structure was subjected to 100 steps of steepest-descent minimization followed by 500 steps of adopted-basis Newton–Raphson minimization. The system was centered at the midpoint of the membrane as the origin and the helical axis of Vpu aligned to the z-axis (the normal to the membrane plane). The system was heated from 0 to 298.15 K. Subsequently, a 60-ns trajectory was generated using Langevin dynamics. Frictional constant of 91 ps−1 was used. The SHAKE (Ryckaert et al. 1977) algorithm was used to constrain bonds involving hydrogen, allowing for a 0.002-ps integration time step.
Free-energy profiles for the tilt angle θ relative to the z-axis were determined by means of umbrella sampling and the weighted histogram analysis method (WHAM). This method has been extensively used to determine free-energy profiles, and as such, the details related to theoretical background and implementation are not included here (see Kumar et al. 1992; Roux 1995). Briefly, umbrella sampling involves applying a bias to restrict sampling along some reaction coordinate. To obtain potential of mean force (PMF) profiles, the reaction coordinate is divided into a number of sampling windows, and in each window, simulations are biased to sample conformations around some value of the reaction coordinate. Statistics are collected in each window, and then WHAM is used to unbias and combine statistics from all windows to generate the desired PMF profiles along the entire range of the reaction coordinate. In this study, the biasing potential was of the following form:\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$ V\left( \theta \right)\, = \,\frac{k}{2}\left( {\theta \, - \,\theta_{0} } \right)^{2} , $$\end{document}where k is the force constant, θ is the tilt angle (defined as the angle subtended by the vector along the helical and the vector normal to the membrane; Fig. 1), and θ0 is the value of the tilt angle around which sampling is desired. Statistics were collected in 36 sampling windows in the range θ0 = [0, 90], with a bin size of 2.5°. In a typical window, a force constant of 250 kcal/mol/deg2 was used. Simulation in each window was 4 ns in length. The first 1 ns was discarded, and the other 3 ns was used to collect statistics.
Protein reorientation in response to hydrophobic mismatch. Schematic representation of the tilting of a TM protein under conditions of hydrophobic mismatch when membrane thickness (D) is smaller than the length of the helix (L). In response to mismatch, the helix tilts relative to the normal of the membrane (here the z direction) such that hydrophobic side chains are buried inside the membrane
Results and Discussion
The variation of the helix tilt as a function of time and membrane thickness, D, is summarized in Fig. 2. For simulations using IMM1 (Fig. 2a) and GBSW (Fig. 2b), the helix quickly tilts in response to positive hydrophobic mismatch, and it thereafter fluctuates about their mean. In the case of IMM1, as the membrane thickness is decreased from 27 to 15 Å, the tilt angle, measured as the time average of the length of the trajectory, increased from 9° to 45° (Fig. 2a). Similarly, in the case of GBSW, as the membrane thickness is decreased from 27 to 15 Å, the tilt angle increased from 18° to 56° (Fig. 2b). For comparison, the corresponding experimentally determined tilt angles are shown in parentheses. These values show very good agreement with mean tilt angles for both IMM1 and GBSW. When cos θ is plotted against D, a roughly linear relationship is revealed, in agreement with Park and Opella (2005) (Fig. 2c). These results corroborate those of other molecular dynamics simulation studies on α-helical TM proteins. Specifically, using an explicit lipid representation for the membrane, Goodyear et al. (2005) observed that for a synthetic peptide initially oriented along the normal of membrane the helix tilts in response to positive hydrophobic mismatch, with an average tilt angle comparable to the experimentally determined value.
Tilt angle time series. Tilt angles are shown as a function of time for Vpu simulated in the IMM1 (a) and GBSW (b) implicit membrane models. Results are shown for simulations where D (membrane thickness) is 15 (black), 19 (blue), 23 (green), and 27 Å (red). Shown on each plot are the average angles. For comparison, experimental tilt angles are shown in parentheses. Shown in (c) is a plot of cos(θ) vs. D determined using results from IMM1 simulation (red), GBSW (black) simulation, and experiment (blue) (Color figure online)
The next phase of the study involved the determination of free-energy profiles as a function of the tilt of the α-helical peptide. By means of umbrella sampling, the free energy as a function of tilt angle was determined for angles ranging from 0° to 90° for membrane thicknesses of 15, 19, 23, and 27 Å, respectively. As mentioned earlier, if we assume that the optimal tilt angles a helix adopts in a membrane represent its energetically favorable orientations according to the free energy computed with the implicit solvent model, then by using these free-energy profiles, we can predict tilt angles. In the case of IMM1, for D = 15, 19, 23, and 27 Å, the absolute minimum was located at 41°, 31°, 17°, and 11°, respectively (Fig. 3a). In the case GBSW, for D = 15, 19, 23, and 27 Å, the absolute minimum was located at 51°, 41°, 29°, and 14°, respectively (Fig. 3b).
Potential of mean force profiles. Free-energy profiles are shown as a function of tilt angle (θ) calculated using the IMM1 (a) and GBSW (b) implicit membrane models. Shown are profiles calculated where D (membrane thickness) is 15 (black), 19 (blue), 23 (green), and 27 Å (red) (Color figure online)
The fact that the results obtained from the IMM1 and GBSW models so closely resemble the experiment is significant because it implies that the implicit models used for this system contained all the necessary physics as it relates to electrostatic and hydrophobic interaction within the membrane region. Although more detailed interactions such as anchoring interactions between peptides and lipids (e.g., tryptophan residues interacting with the polar head group region; de Planque et al. 2003) can compete with the fundamental free-energy driving force captured in IMM1 and GBSW, we believe that our study shows that IMM1 and GBSW, at least for systems as simple as the one used in this project, certainly embody the necessary functional form and parameters to describe the tilt process as a function of positive hydrophobic mismatch. IMM1 and GBSW thus prove promising for membrane proteins, given the savings in computer time, while at the same time accurately recreating the solvent effects of the cell membrane and the surrounding aqueous regions.
A comparison between IMM1 and GBSW reveals that they yield qualitatively similar results. They were both able to recreate the tilting of the α-helical TM peptide Vpu as function of hydrophobic mismatch that was observed experimentally, and they both exhibit a linear relationship between cos θ and D, in agreement with Park and Opella (2005). Quantitatively, however, IMM1 and GBSW give slightly different tilt angles for a given membrane thickness. In general, IMM1 gave values lower than the experimental observed values, while GBSW gave higher values (Fig 2). This is also reflected in the free-energy profile, where qualitatively the models mirrored each other, but quantitatively there were slight difference (GBSW predicted the minima to exist at higher θ for any given D; Fig. 3). The PMFs for IMM1 and GBSW also differ as the tilt angle approaches 80°. For IMM1, the barrier to tilting beyond the minima is larger as the membrane thickness decreases. The trend is reversed slightly in the case of GBSW, but in this region, the difference amounts to a few kcal/mol.
In terms of the performance of the algorithms, simulations using IMM1 were on average four times as fast as GBSW. This can be attributed to the fact that GBSW involves the time-consuming calculation of the accessible surface area, which is used to approximate the nonpolar contribution to the solvation energy.
An intriguing knowledge-based potential by Ulmschneider et al. (2005) was derived and recently applied to the Vpu helix (Ulmschneider et al. 2006). The basic hypothesis of that study was that each residue on its own would prefer a certain depth inside the membrane bilayer. Potentials of mean force along the membrane normal were derived for each amino acid by fitting Gaussian functions to residue position distributions from known TM structures. The individual potentials agreed well with experimental and theoretical considerations. The resulting implicit membrane potential was tested on various membrane proteins as well as single trans-membrane helices. All membrane proteins were found to be at an energy minimum when correctly inserted into the membrane. The results of that study qualitatively match those obtained in our study.
Conclusions
In our computational test, the use of two widely used implicit membrane models, IMM1 and GBSW, was sufficient to reproduce the experimentally measured tilting (Park and Opella 2005) of an α-helical peptide that occurs in response to positive hydrophobic mismatch. Umbrella sampling-based free-energy profile calculations to obtain the distribution of the peptide tilt angles relative to the membrane proved to be a reliable means of gaining information concerning the free energy of the system.
ATF was supported by the National Science Foundation Graduate Fellowship program. IA acknowledges support from the National Science Foundation Career award (CHE-0918817).
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
ReferenceBartelsCKarplusMMultidimensional adaptive umbrella sampling: applications to main chain and side chain peptide conformationsJ Comput Chem1997181450146210.1002/(SICI)1096-987X(199709)18:12<1450::AID-JCC3>3.0.CO;2-IBennettCHEfficient estimation of free-energy differences from Monte-Carlo dataJ Comput Phys19762224526810.1016/0021-9991(76)90078-4BranniganGLinLCLBrownFLHImplicit solvent simulation models for biomembranesEur Biophys J Biophys Lett200635104124PlanqueMRRBonevBBDemmersJAAGreathouseDVKoeppeREIISeparovicFWattsAKillianJAInterfacial anchor properties of tryptophan residues in transmembrane peptides can dominate over hydrophobic matching effects in peptide–lipid interactionsBiochemistry2003425341534810.1021/bi027000r12731875FeigMBrooksCLRecent advances in the development and application of implicit solvent models in biomolecule simulationsCurr Opin Struct Biol20041421722410.1016/j.sbi.2004.03.00915093837ForrestLRSansomMSPMembrane simulations: bigger and betterCurr Opin Struct Biol20001017418110.1016/S0959-440X(00)00066-X10753807GoodyearDJSharpeSGrantCWMMorrowMRMolecular dynamics simulation of transmembrane polypeptide orientational fluctuationsBiophys J20058810511710.1529/biophysj.104.04750615489306ImWBrooksCLInterfacial folding and membrane insertion of designed peptides studied by molecular dynamics simulationsProc Natl Acad Sci U S A20051026771677610.1073/pnas.040813510215860587ImWFeigMBrooksCLAn implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteinsBiophys J2003852900291810.1016/S0006-3495(03)74712-214581194JarzynskiCNonequilibrium equality for free energy differencesPhys Rev Lett1997782690269310.1103/PhysRevLett.78.2690KarplusMMcCammonJAMolecular dynamics simulations of biomoleculesNat Struct Biol2002964665210.1038/nsb0902-64612198485KillianJAHydrophobic mismatch between proteins and lipids in membranesBiochim Biophys Acta199813764014169805000KirkwoodJGStatistical mechanics of fluid mixturesJ Chem Phys1935330010.1063/1.1749657KumarSBouzidaDSwendsenRHKollmanPARosenbergJMThe weighted histogram analysis method for free-energy calculations on biomolecules. I. The methodJ Comput Chem1992131011102110.1002/jcc.540130812LazaridisTEffective energy function for proteins in lipid membranesProteins20035217619210.1002/prot.1041012833542LazaridisTKarplusMEffective energy function for proteins in solutionProteins19993513315210.1002/(SICI)1097-0134(19990501)35:2<133::AID-PROT1>3.0.CO;2-N10223287LevyRMGallicchioEComputer simulations with explicit solvent: recent progress in the thermodynamic decomposition of free energies and in modeling electrostatic effectsAnnu Rev Phys Chem19984953156710.1146/annurev.physchem.49.1.5319933909MottamalMZhangJLazaridisTEnergetics of the native and non-native states of the glycophorin transmembrane helix dimerProteins200662996100910.1002/prot.2084416395713ParkSHOpellaSJTilt angle of a trans-membrane helix is determined by hydrophobic mismatchJ Mol Biol200535031031810.1016/j.jmb.2005.05.00415936031ParkSHMrseAANevzorovAAMeslehMFOblatt-MontalMMontalMOpellaSJThree-dimensional structure of the channel-forming trans-membrane domain of virus protein “u” (Vpu) from HIV-1J Mol Biol200333340942410.1016/j.jmb.2003.08.04814529626PastorRWMolecular dynamics and Monte Carlo simulations of lipid bilayersCurr Opin Struct Biol1994448649210.1016/S0959-440X(94)90209-7RouxBCalculation of the potential of mean force using computer simulationsComput Phys Commun19959127528210.1016/0010-4655(95)00053-IRouxBSimonsonTImplicit solvent modelsBiophys Chem19997812010.1016/S0301-4622(98)00226-917030302RyckaertJPCiccottiGBerendsenHJCNumerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanesJ Comput Phys19772332734110.1016/0021-9991(77)90098-5TorrieGMValleauJPNon-physical sampling distributions in Monte-Carlo free-energy estimation—umbrella samplingJ Comput Phys19772318719910.1016/0021-9991(77)90121-8UlmschneiderMBSansomMSPNolaAProperties of integral membrane protein structures: derivation of an implicit membrane potentialProteins20055925226510.1002/prot.2033415723347UlmschneiderMBSansomMSPNolaAEvaluating tilt angles of membrane-associated helices: comparison of computational and NMR techniquesBiophys J2006901650166010.1529/biophysj.105.06536716339877WhiteSHHeijneGTransmembrane helices before, during, and after insertionCurr Opin Struct Biol20051537838610.1016/j.sbi.2005.07.00416043344ZwanzigRWHigh-temperature equation of state by a perturbation method. I. Nonpolar gasesJ Chem Phys195422142010.1063/1.1740193