Dielectric Continuum Models for Hydration Effects on Peptide Conformational Transitions

Quantifying the Molecular Polarization Response of Liquid Water Interfaces at Heterogeneously Charged Surfaces

The hydration shells of proteins mediate interactions, such as small molecule binding, that are vital to their biological function or in some cases their dysfunction. However, even when the structure of a protein is known, the properties of its hydration environment cannot be easily predicted due to the complex interplay between protein surface heterogeneity and the collective structure of water's hydrogen bonding network. This manuscript presents a theoretical study of the influence of surface charge heterogeneity on the polarization response of the liquid water interface. We focus our attention on classical point charge models of water, where the polarization response is limited to molecular reorientation. We introduce a new computational method for analyzing simulation data that is capable of quantifying water's collective polarization response and determining the effective surface charge distribution of hydrated surfaces over atomistic length scales. To illustrate the utility of this method, we present the results of molecular dynamics simulations of liquid water in contact with a heterogeneous model surface and the CheY protein.

Biomolecular electrostatics and solvation: a computational perspective

An understanding of molecular interactions is essential for insight into biological systems at the molecular scale. Among the various components of molecular interactions, electrostatics are of special importance because of their long-range nature and their influence on polar or charged molecules, including water, aqueous ions, proteins, nucleic acids, carbohydrates, and membrane lipids. In particular, robust models of electrostatic interactions are essential for understanding the solvation properties of biomolecules and the effects of solvation upon biomolecular folding, binding, enzyme catalysis, and dynamics. Electrostatics, therefore, are of central importance to understanding biomolecular structure and modeling interactions within and among biological molecules. This review discusses the solvation of biomolecules with a computational biophysics view toward describing the phenomenon. While our main focus lies on the computational aspect of the models, we provide an overview of the basic elements of biomolecular solvation (e.g. solvent structure, polarization, ion binding, and non-polar behavior) in order to provide a background to understand the different types of solvation models.

Computational methods for biomolecular electrostatics

An understanding of intermolecular interactions is essential for insight into how cells develop, operate, communicate, and control their activities. Such interactions include several components: contributions from linear, angular, and torsional forces in covalent bonds, van der waals forces, as well as electrostatics. Among the various components of molecular interactions, electrostatics are of special importance because of their long range and their influence on polar or charged molecules, including water, aqueous ions, and amino or nucleic acids, which are some of the primary components of living systems. Electrostatics, therefore, play important roles in determining the structure, motion, and function of a wide range of biological molecules. This chapter presents a brief overview of electrostatic interactions in cellular systems, with a particular focus on how computational tools can be used to investigate these types of interactions.

Finding the transition state without initial guess: The growing string method for Newton trajectory to isomerization and enantiomerization reaction of alanine dipeptide and poly(15)alanine

We report a new, high-dimensional application of a method for finding a transition state (TS) between a reactant and a product on the potential energy surface: the search of a growing string along a reaction path defined by any Newton trajectory in combination with the Berny method (Quapp, J Chem Phys (2005), 122, 174106; we have provided this algorithm on a web page). Two given minima are connected by a one-dimensional, but usually curvilinear reaction coordinate. It leads to the TS region. The application of the method to alanine dipeptide finds the TS of the isomerisation C(7 ax) --> C(5), some TSs of the enantiomerisation of C(7 ax) from L-form to quasi-D-form, and it finds the TS region of a transition of a partly unfolded, bent structure which turns back into a mainly alpha-helix in the Ac(Ala)(15)NHMe polyalanine (all at the quantum mechanical level B3LYP/6-31G: the growing string calculation is interfaced with the Gaussian03 package). The formation or dissolvation of some alpha- or 3(10)-hydrogen bonds of the helix are discussed along the TS pathway, as well as the case of an enantiomer at the central residue of the helix.

Structural details, pathways, and energetics of unfolding apomyoglobin

Protein folding is often difficult to characterize experimentally because of the transience of intermediate states, and the complexity of the protein-solvent system. Atomistic simulations, which could provide more detailed information, have had to employ highly simplified models or high temperatures, to cope with the long time scales of unfolding; direct simulation of folding is even more problematic. We report a fully atomistic simulation of the acid-induced unfolding of apomyoglobin in which the protonation of acidic side-chains to simulate low pH is sufficient to induce unfolding at room temperature with no added biasing forces or other unusual conditions; and the trajectory is validated by comparison to experimental characterization of intermediate states. Novel insights provided by their analysis include: characterization of a dry swollen globule state forming a barrier to initial unfolding or final folding; observation of cooperativity in secondary and tertiary structure formation and its explanation in terms of dielectric environments; and structural details of the intermediate and the completely unfolded states. These insights involve time scales and levels of structural detail that are presently beyond the range of experiment, but come within reach through the simulation methods described here. An implicit solvation model is used to analyze the energetics of protein folding at various pH and ionic strength values, and a reasonable estimate of folding free energy is obtained. Electrostatic interactions are found to disfavor folding.

Ion permeation through the gramicidin channel: atomically detailed modeling by the Stochastic Difference Equation

Atomically detailed descriptions of ionic solution, membrane, and the gramicidin channel are used to compute molecular dynamics trajectories of ion permeation. The microsecond trajectories are calculated with the Stochastic Difference Equation (SDE), which provides approximate solutions to the equations of motions (with filtered high-frequency modes) of extended timescales. The relative permeations of lithium, sodium, and potassium are estimated by using a novel, kinetic cycle protocol and are compared with experiment. The transport through native gramicidin and one fluoro-valine variant is considered as well. Qualitative agreement between theory and experiment is obtained. The faster permeation rate of sodium compared to lithium is reproduced in the calculations. The calculations also reproduce the slower diffusion through a gramicidin with fluorinated valine compared to native gramicidin. The calculations are inconclusive about the relative rates of potassium and sodium. The experiment suggests that potassium permeates more quickly. We directly probe the kinetics of a biophysical process at a relevant time window without reducing the atomically detailed description of the system. The calculations were able to capture subtle balances between binding and diffusion that determine permeation rates. The same model gave the correct ordering of diffusion rates for cases in which electrostatic binding has opposite effects and must be supplemented by dynamic factors. Diffusion rates are faster when favorable electrostatic interactions of ions in the channel (compared to the solvent) are observed. Studies of a gramicidin variant suggest an opposite effect, in which permeation is faster for the less polar channel, indicating dynamic effects. Although both trends can be explained qualitatively, it is not possible to predict (before doing the SDE calculations) which factor is more important.

Accelerated Poisson-Boltzmann calculations for static and dynamic systems

We report here an efficient implementation of the finite difference Poisson-Boltzmann solvent model based on the Modified Incomplete Cholsky Conjugate Gradient algorithm, which gives rather impressive performance for both static and dynamic systems. This is achieved by implementing the algorithm with Eisenstat's two optimizations, utilizing the electrostatic update in simulations, and applying prudent approximations, including: relaxing the convergence criterion, not updating Poisson-Boltzmann-related forces every step, and using electrostatic focusing. It is also possible to markedly accelerate the supporting routines that are used to set up the calculations and to obtain energies and forces. The resulting finite difference Poisson-Boltzmann method delivers efficiency comparable to the distance-dependent dielectric model for a system tested, HIV Protease, making it a strong candidate for solution-phase molecular dynamics simulations. Further, the finite difference method includes all intrasolute electrostatic interactions, whereas the distance dependent dielectric calculations use a 15-A cutoff. The speed of our numerical finite difference method is comparable to that of the pair-wise Generalized Born approximation to the Poisson-Boltzmann method.

Distinguishing native conformations of proteins from decoys with an effective free energy estimator based on the OPLS all-atom force field and the Surface Generalized Born solvent model

Protein decoy data sets provide a benchmark for testing scoring functions designed for fold recognition and protein homology modeling problems. It is commonly believed that statistical potentials based on reduced atomic models are better able to discriminate native-like from misfolded decoys than scoring functions based on more detailed molecular mechanics models. Recent benchmark tests on small data sets, however, suggest otherwise. In this work, we report the results of extensive decoy detection tests using an effective free energy function based on the OPLS all-atom (OPLS-AA) force field and the Surface Generalized Born (SGB) model for the solvent electrostatic effects. The OPLS-AA/SGB effective free energy is used as a scoring function to detect native protein folds among a total of 48,832 decoys for 32 different proteins from Park and Levitt's 4-state-reduced, Levitt's local-minima, Baker's ROSETTA all-atom, and Skolnick's decoy sets. Solvent electrostatic effects are included through the Surface Generalized Born (SGB) model. All structures are locally minimized without restraints. From an analysis of the individual energy components of the OPLS-AA/SGB energy function for the native and the best-ranked decoy, it is determined that a balance of the terms of the potential is responsible for the minimized energies that most successfully distinguish the native from the misfolded conformations. Different combinations of individual energy terms provide less discrimination than the total energy. The results are consistent with observations that all-atom molecular potentials coupled with intermediate level solvent dielectric models are competitive with knowledge-based potentials for decoy detection and protein modeling problems such as fold recognition and homology modeling.

A classical and ab initio study of the interaction of the myosin triphosphate binding domain with ATP

We used classical molecular mechanics (MM) simulations and quantum mechanical (QM) structural relaxations to examine the active site of myosin when bound to ATP. Two conformations of myosin have been determined by x-ray crystallography. In one, there is no direct interaction between switch 2 and the nucleotide (open state). In the other (closed state), the universally conserved switch 2 glycine forms a hydrogen bond with a gamma-phosphate oxygen. MM simulations indicate that the two states are thermodynamically stable and allow us to investigate the extent to which the P-loop, switch 1, and switch 2 are involved in hydrolysis. We find that the open structure has a higher affinity for ATP than the closed structure, and that ATP is distorted toward a transition state by interactions with the protein. We also examine how the structure of the binding site changes with either MgATP or CaATP as the nucleotide in myosin in the open conformer. Our analyses suggest that higher CaATPase rates occur because the leaving phosphate (P(i)) group is more weakly bound and dissociation occurs faster. Finally, we validate the use of a particular formulation of a QM methodology (Car-Parrinello) to further refine the structures of the active site.

Artificial protein cavities as specific ligand-binding templates: characterization of an engineered heterocyclic cation-binding site that preserves the evolved specificity of the parent protein

Cavity complementation has been observed in many proteins, where an appropriate small molecule binds to a cavity-forming mutant. Here, the binding of compounds to the W191G cavity mutant of cytochrome c peroxidase is characterized by X-ray crystallography and binding thermodynamics. Unlike cavities created by removal of hydrophobic side-chains, the W191G cavity does not bind neutral or hydrophobic compounds, but displays a strong specificity for heterocyclic cations, consistent with the role of the protein to stabilize a tryptophan radical at this site. Ligand dissociation constants for the protonated cationic state ranged from 6 microM for 2-amino-5-methylthiazole to 1 mM for neutral ligands, and binding was associated with a large enthalpy-entropy compensation. X-ray structures show that each of 18 compounds with binding behavior bind specifically within the artificial cavity and not elsewhere in the protein. The compounds make multiple hydrogen bonds to the cavity walls using a subset of the interactions seen between the protein and solvent in the absence of ligand. For all ligands, every atom that is capable of making a hydrogen bond does so with either protein or solvent. The most often seen interaction is to Asp235, and most compounds bind with a specific orientation that is defined by their ability to interact with this residue. Four of the ligands do not have conventional hydrogen bonding atoms, but were nevertheless observed to orient their most polar CH bond towards Asp235. Two of the larger ligands induce disorder in a surface loop between Pro190 and Asn195 that has been identified as a mobile gate to cavity access. Despite the predominance of hydrogen bonding and electrostatic interactions, the small variation in observed binding free energies were not correlated readily with the strength, type or number of hydrogen bonds or with calculated electrostatic energies alone. Thus, as with naturally occurring binding sites, affinities to W191G are likely to be due to a subtle balance of polar, non-polar, and solvation terms. These studies demonstrate how cavity complementation and judicious choice of site can be used to produce a protein template with an unusual ligand-binding specificity.

Continuum solvent molecular dynamics study of flexibility in interleukin-8

Generalized Born continuum solvent methods have been shown to provide a reasonable description of the equilibrium thermodynamics of aqueous solvation in a variety of applications to peptides, proteins, and nucleic acids. Here we study the performance of these methods in molecular dynamics simulations of interleukin-8, comparing nanosecond-length explicit solvent simulations with those using the generalized Born model. In general, the simulations show similar results, although movement away from the initial NMR-determined structure and average fluctuations about the mean are slightly higher for the continuum solvent results. In both simulations, the two helices that are packed on top of the core sheet move closer together, resulting in a structure that more closely resembles the X-ray structure. Principal-component (quasiharmonic) analysis is used to analyze the motions of these helices in both of the simulations and in the NMR ensemble of structures. Prospects for making more general use of continuum solvent models in protein dynamics simulations are discussed.

Applications of neural network prediction of conformational states for small peptides from spectra and of fold classes

Electronic structures of small peptides were calculated 'ab initio' with the help of Density Functional Theory (DFT) and molecular dynamics that rendered a set of conformational states of the peptides. For the structures of these states it was possible to derive atomic polar tensors that allowed us to construct vibrational spectra for each of the conformational states with low energy. From the spectra, neural networks could be trained to distinguish between the various states and thus be able to generate a larger set of relevant structures and their relation to secondary structures of the peptides. The calculations were done both with solvent atoms (up to ten water molecules) and without, and hence the neural networks could be used to monitor the influence of the solvent on hydrogen bond formation. The calculations at this stage only involved very short peptide fragments of a few alanine amino acids but already at this stage they could be compared with reasonable agreements to experiments. The neural networks are shown to be good in distinguishing the different conformers of the small alanine peptides, especially when in the gas phase. Also the task of predicting protein fold-classes, defined from line-geometry, seems promising.

Neural-network analysis of the vibrational spectra of N-acetyl L-alanyl N'-methyl amide conformational states

Density-functional theory (DFT) calculations utilizing the Becke 3LYP hybrid functional have been carried out for N-acetyl L-alanine N'-methylamide and examined with respect to the effect of water on the structure, the vibrational frequencies, vibrational absorption (VA), vibrational circular dichroism (VCD), Raman spectra, and Raman optical activity (ROA) intensities. The large changes due to hydration in the structures, and the relative stability of the conformer, reflected in the VA, VCD, Raman spectra, and ROA spectra observed experimentally, are reproduced by the DFT calculations. A neural network has been constructed for reproducing the inverse scattering data (we infer the structural coordinates from spectroscopic data) that the DFT method could produce. The purpose of the network has also been to generate the large set of conformational states associated with each set of spectroscopic data for a given conformer of the molecule by interpolation. Finally the neural network performances are used to monitor a sensitivity analysis of the importance of secondary structures and the influence of the solvent. The neural network is shown to be good in distinguishing the different conformers of the small alanine peptide, especially in the gas phase.

Solvation energy density occlusion approximation for evaluation of desolvation penalties in biomolecular interactions

In calculations involving many displacements of an interacting pair of biomolecules, such as brownian dynamics, the docking of a substrate/ligand to an enzyme/receptor, or the screening of a large number of ligands as prospective inhibitors for a particular receptor site, there is a need for rapid evaluation of the desolvation penalties of the interacting pair. Although continuum electrostatic treatments with distinct dielectric constants for solute and solvent provide an account of the electrostatics of solvation and desolvation, it is necessary to re-solve the Poisson equation, at considerable computational cost, for each displacement of the interacting pair. We present a new method that uses a formulation of continuum electrostatic solvation in terms of the solvation energy density and approximates desolvation in terms of the occlusion of this density. We call it the SEDO approximation. It avoids the need to re-solve the Poisson equation, as desolvation is now estimated by an integral over the occluded volume. Test calculations are presented for some simple model systems and for some real systems that have previously been studied using the Poisson equation approach: MHC class I protein-peptide complexes and a congeneric series of human immunodeficiency virus type 1 (HIV-1) protease--ligand complexes. For most of the systems considered, the trends and magnitudes of the desolvation component of interaction energies obtained using the SEDO approximation are in reasonable correlation with those obtained by re-solving the Poisson equation. In most cases, the error introduced by the SEDO approximation is much less than that of the often-used test-charge approximation for the charge-charge components of intermolecular interactions. Proteins 2001;43:12-27.

Molecular dynamics study of the energetic, mechanistic, and structural implications of a closed phosphate tube in ncd

The switch 1 region of myosin forms a lid over the nucleotide phosphates as part of a structure known as the phosphate-tube. The homologous region in kinesin-family motors is more open, not interacting with the nucleotide. We used molecular dynamics (MD) simulations to examine a possible displacement of switch 1 of the microtubule motor, ncd, from the open conformation to the closed conformation seen in myosin. MD simulations were done of both the open and the closed conformations, with either MgADP or MgATP at the active site. All MD structures were stable at 300 K for 500 ps, implying that the open and closed conformers all represented local minima on a global free energy surface. Free energy calculations indicated that the open structure was energetically favored with MgADP at the active site, suggesting why only the open structure has been captured in crystallographic work. With MgATP, the closed and open structures had roughly equal energies. Simulated annealing MD showed the transformation from the closed phosphate-tube ncd structure to an open configuration. The MD simulations also showed that the coordination of switch 1 to the nucleotide dramatically affected the position of both the bound nucleotide and switch 2 and that a closed phosphate-tube may be necessary for catalysis.

Generalized born models of macromolecular solvation effects

It would often be useful in computer simulations to use a simple description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation, and approximations to such models that avoid the need to solve the Poisson equation are attractive because of their computational efficiency. Here we give an overview of one such approximation, the generalized Born model, which is simple and fast enough to be used for molecular dynamics simulations of proteins and nucleic acids. We discuss its strengths and weaknesses, both for its fidelity to the underlying continuum model and for its ability to replace explicit consideration of solvent molecules in macromolecular simulations. We focus particularly on versions of the generalized Born model that have a pair-wise analytical form, and therefore fit most naturally into conventional molecular mechanics calculations.

Effective energy function for proteins in solution

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.

Nucleic acid base-pairing and N-methylacetamide self-association in chloroform: affinity and conformation

A recently developed computational method, 'mining minima', is used to examine the hydrogen-bonding interactions of nucleic acid base-pairs and of the N-methylacetamide homodimer in chloroform. The mining minima algorithm aggressively samples molecular conformations, identifies the most important local minima, and computes their contributions to the overall free energy of the system. Here, the CHARMM 98 parameter set is used for the potential energy and the generalized Born/surface area solvent model is used to account for the influence of the solvent. Good agreement with experiment is obtained for the non-covalent binding affinities of a series of complexes. The computational approach used here is applicable to a range of molecular systems.

Continuum solvent studies of the stability of RNA hairpin loops and helices

We apply continuum solvent models to investigate the relative stability of various conformational forms for two RNA sequences, GGAC(UUCG)GUCC and GGUG(UGAA)CACC. In the first part, we compare alternate hairpin conformations to explore the reliability of these models to discriminate between different local conformations. A second part looks at the hairpin-duplex conversion for the UUCG sequence, identifying major contributors to the thermodynamics of a much large scale transition. Structures were taken as snapshots from multi-nanosecond molecular dynamics simulations computed in a consistent fashion using explicit solvent and with long-range electrostatics accounted for using the Particle-Mesh Ewald procedure. The electrostatic contribution to solvation energies were computed using both a finite-difference Poisson-Boltzmann (PB) model and a pairwise Generalized Born model; non-electrostatic contributions were estimated with a surface-area dependent term. To these solvation free energies were added the mean solute internal energies (determined from a molecular mechanics potential) and estimates of the solute entropy (from a harmonic analysis). Consistent with experiment and with earlier solvated molecular dynamics simulations, the UUCG hairpin was found to prefer conformers close to a recent NMR structure determination in preference to those from an earlier NMR study. Similarly, results for the UGAA hairpin favored an NMR-derived structure over that to be expected for a generic GNRA hairpin loop. Experimental free energies are not known for the hairpin/duplex conversion, but must be close to zero since hairpins are seen in solution and duplexes in crystals; out calculations find a value near zero and illustrate the expected interplay of solvation, salt effects and entropy in affecting this equilibrium.

Molecular dynamics and continuum solvent studies of the stability of polyG-polyC and polyA-polyT DNA duplexes in solution

Molecular dynamics simulation in explicit solvent and continuum solvent models are applied to investigate the relative stability of A- and B-form helices for two DNA sequences, dA10-dT10 and dG10-dC10 in three structural forms. One structural form is based on an unrestrained molecular dynamics (MD) trajectory starting from a canonical B-DNA structure, the second is based on a MD trajectory starting in a canonical B-DNA structure with the sugars constrained to be C2'-endo and the third simulation started from a canonical A-DNA structure with the sugars constrained to C3'-endo puckers. For the energetic analysis, structures were taken as snapshots from nanosecond length molecular dynamics simulations computed in a consistent fashion in explicit solvent, applying the particle mesh Ewald method and the Cornell et al. force field. The electrostatic contributions to solvation free energies are computed using both a finite-difference Poisson-Boltzmann model and a pairwise Generalized Born model. The non-electrostatic contributions to the solvation free energies are estimated with a solvent accessible surface area dependent term. To estimate the gas phase component of the relative free energy between the various structures, the mean solute internal energies (determined with the Cornell et al. molecular mechanics potential including all pairwise interactions within the solute) and estimates of the solute entropy (using a harmonic approximation) were used. Consistent with experiment, the polyG-polyC (GC) structures are found to be much more A-phillic than the polyA-polyT (AT) structures, the latter being quite A-phobic. The dominant energy components responsible for this difference comes from the internal and van der Waal energies. A perhaps less appreciated difference between the GC and AT rich sequences is suggested by the calculated salt dependence which demonstrates a significantly enhanced ability to drive GC rich sequences towards an A-form structure compared to AT rich sequences. In addition to being A-phobic, the AT structure also has a noticably larger helical repeat than GC and other mixed sequence duplexes, consistent with experiment. Analysis of the average solvent density from the trajectories shows hydration patterns in qualitative agreement with experiment and previous theoretical treatments.

Hydration and conformational equilibria of simple hydrophobic and amphiphilic solutes

We consider whether the continuum model of hydration optimized to reproduce vacuum-to-water transfer free energies simultaneously describes the hydration free energy contributions to conformational equilibria of the same solutes in water. To this end, transfer and conformational free energies of idealized hydrophobic and amphiphilic solutes in water are calculated from explicit water simulations and compared to continuum model predictions. As benchmark hydrophobic solutes, we examine the hydration of linear alkanes from methane through hexane. Amphiphilic solutes were created by adding a charge of +/-1e to a terminal methyl group of butane. We find that phenomenological continuum parameters fit to transfer free energies are significantly different from those fit to conformational free energies of our model solutes. This difference is attributed to continuum model parameters that depend on solute conformation in water, and leads to effective values for the free energy/surface area coefficient and Born radii that best describe conformational equilibrium. In light of these results, we believe that continuum models of hydration optimized to fit transfer free energies do not accurately capture the balance between hydrophobic and electrostatic contributions that determines the solute conformational state in aqueous solution.

Thermodynamics of a reverse turn motif. Solvent effects and side-chain packing

The linear pentapeptide, Ala-Tyr-cis-Pro-Tyr-Asp-NMA (AYPYD) is known to have a significant population of type VI turn conformers in aqueous solvent. We have carried out theoretical studies of the conformational energetics of this peptide using a potential of mean force (PMF) consisting of the AMBER/OPLS empirical potential energy function, a macroscopic electrostatic model of polar solvation, and a surface area-based model of non-polar solvation. Conformers were taken from molecular dynamics simulations reported elsewhere, or generated by a random search method reported here. The chain entropy of folding was calculated by a systematic search of accessible dihedral angle space. The intra-peptide component was found to strongly favor folding and was nearly cancelled by the polar solvation term which disfavored folding. The non-polar solvation term had little effect. Fluctuations about the average value of the PMF were small and in accord with estimates from a simple harmonic model. When applied to conformers generated by a random search, the PMF selected a conformer close to the NMR-determined structure as the lowest energy conformer. The conformer with the second-lowest energy was extended, but was found to fold rapidly to the turn state in a subsequent molecular dynamics study, and may be an important state on the folding-unfolding pathway. Averages of the PMF were combined with the entropy estimates to provide an estimate of the free energy of folding that is in reasonable agreement with experimental results. In terms of the interplay between backbone electrostatic interactions and the packing of apolar side-chains, this peptide provides a model for the energetics of protein folding, and therefore makes a useful test case for calculations.