| MC-PHS: A Monte Carlo Implementation of the Primary Hydration Shell for Protein Folding and Design Biophysical Journal, Volume 84, Issue 2, 1 February 2003, Pages 805-815 Alex Kentsis, Mihaly Mezei and Roman Osman Abstract A primary hydration shell (PHS) approach is developed for Monte Carlo simulations of conformationally rich macromolecular systems in an environment that efficiently captures principal solvation effects. It has been previously demonstrated that molecular dynamics using PHS is an efficient method to study peptide structure and dynamics in aqueous solution. Here, we extend the PHS approach to Monte Carlo simulations, whereby a stable shell of water molecules is maintained with a flexible, nonspherical, half-harmonic potential, tuned to maintain a constant restraining energy, with the difference between the restraint and shell energies used to dynamically adjust the shell radius. Examination of the shell and system size dependence of the restraining potential reveals its robustness. Moreover, its suitability for biomolecular simulations is evaluated using small spheres of water, hydration properties of small biological molecules, and configurational sampling of -hairpin pentapeptide YPGDV. This method, termed MC-PHS, appears to provide efficient representation of dominant solvation effects and should prove useful in the study of protein folding and design. Abstract | Full Text | PDF (537 kb) |
| Application of the Primary Hydration Shell Approach to Locally Enhanced Sampling Simulated Annealing: Computer Simulation of Thyrotropin-Releasing Hormone in Water Biophysical Journal, Volume 79, Issue 1, 1 July 2000, Pages 66-79 Avia Rosenhouse-Dantsker and Roman Osman Abstract A unified model of simulated annealing with locally enhanced sampling (LES) in a primary hydration shell (PHS) aqueous environment is developed and tested by predicting the structure of the tripeptide thyrotropin-releasing hormone (TRH) in solution. The model extends the formulation of the restraining force in the PHS method as a function of temperature, number of copies in the LES method, and shell thickness. The dependence of the restraining force on temperature can be shown to follow the relationship , which can be derived from the expression for kinetic energy in molecular dynamics simulations. The calibration of the restraining force for different simulation conditions reveals the dependence of and on the number of copies in the LES method and the thickness of the PHS. The predicted structure of TRH is in very good agreement with results from NMR experiments and from a 10-ns PHS simulation at 300K. The method promises to be useful in predicting structure of peptides and proteins in an aqueous environment. Abstract | Full Text | PDF (1168 kb) |
| High Apparent Dielectric Constants in the Interior of a Protein Reflect Water Penetration Biophysical Journal, Volume 79, Issue 3, 1 September 2000, Pages 1610-1620 John J. Dwyer, Apostolos G. Gittis, Daniel A. Karp, Eaton E. Lattman, Daniel S. Spencer, Wesley E. Stites and Bertrand García-Moreno E. Abstract A glutamic acid was buried in the hydrophobic core of staphylococcal nuclease by replacement of Val-66. Its pK was measured with equilibrium thermodynamic methods. It was 4.3 units higher than the pK of Glu in water. This increase was comparable to the ΔpK of 4.9 units measured previously for a lysine buried at the same location. According to the Born formalism these ΔpK are energetically equivalent to the transfer of a charged group from water to a medium of dielectric constant of 12. In contrast, the static dielectric constants of dry protein powders range from 2 to 4. In the crystallographic structure of the V66E mutant, a chain of water molecules was seen that hydrates the buried Glu-66 and links it with bulk solvent. The buried water molecules have never previously been detected in >20 structures of nuclease. The structure and the measured energetics constitute compelling and unprecedented experimental evidence that solvent penetration can contribute significantly to the high apparent polarizability inside proteins. To improve structure-based calculations of electrostatic effects with continuum methods, it will be necessary to learn to account quantitatively for the contributions by solvent penetration to dielectric effects in the protein interior. Abstract | Full Text | PDF (379 kb) |
Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 7, L49-L51, 1 April 2007
doi:10.1529/biophysj.106.103010
Biophysical Letters
Mehdi Bagheri Hamaneh and Matthias Buck
, 
Address reprint requests and inquiries to Matthias Buck.Water plays a crucial role for the stability, dynamics, and function of proteins 1,2. For this reason molecular dynamics (MD) simulations must account for the effects that this solvent has, both on protein structure and on protein dynamics. Using a box full of explicitly represented water molecules with periodic boundary conditions (PBC) is the most common way to achieve this goal. However, even in favorable cases the majority of the calculations in such simulations involve water molecules alone. This high computational cost has resulted in several alternative approaches to account for water and water-protein interactions. They can be broadly classified into three categories. In one class, water molecules are entirely replaced by an implicit solvent model (e.g., Lazaridis and Karplus 3 and Im et al. 4). The second category comprises methods that use only a thin shell of explicit waters around the protein 5,6,7. A hybrid explicit/implicit method (a thin shell of explicit waters surrounded by an implicit solvent representation) is used in the third class of models 8,9.
The “primary hydration shell” (PHS) method, developed by Beglov and Roux 5 and implemented in CHARMM, is a good representative of the second category of solvation approaches. In this method, the protein is solvated by a two- to three-layer thick shell of explicit waters, and a half-harmonic constraint is applied to the water molecules if their distance from the nearest protein atom is larger than a certain value (see Supplementary Material for details). Importantly, the force is applied toward the closest protein atom, which results in a water-shell that has the same shape as the protein and can adapt to protein conformational changes. Although the PHS method has been used, for example, for simulated annealing of a tripeptide conformation 10, reports of its application have been few. An explanation for this lack of popularity is that thin water shell models have been thought to lead to unrealistic water and protein behavior. Here we show that the PHS approach leads to acceptable water and especially protein behavior, and thus deserves more attention.
In this study we compare the PHS approach 5 with a simulation using the classical method of solvating the protein in a box full of explicitly represented waters under periodic boundary conditions (PBC). Comparisons are also made with results from an implicit solvent approach, using a recently developed generalized Born method with a simple smoothing function (GBSW) 4. Finally, we present an approach involving a simple harmonic constraint by use of the GEO facility in CHARMM. This approach is similar to the primary hydration shell method but saves on computer time as distances and forces are not calculated with respect to the nearest protein atom but relative to three perpendicular principal axes that follow the protein frame. Such a treatment of forces also allows waters to follow possible changes in protein shape (see Supplementary Material ).
The simulations were carried out on hen lysozyme (Protein Data Bank identifier 6LYT), an ellipsoid-shaped protein that has been subject to extensive experimental 11 and computational 12 studies. The simulation details are described in the Supplementary Material and are briefly summarized as follows: For all simulations involving explicit waters, lysozyme was first immersed in a cubic box of water with a side length of 61.5Å. In the reference simulation PBCs were used together with the particle-mesh Ewald method. Eight chloride atoms neutralized the system’s net charge. All solvent molecules with distances from protein atoms <2.8Å were deleted, eliminating solvent-protein overlap and resulting in 5749 waters. In the case of the PHS approach and its variant employing a simpler harmonic restraint (GEO, described in Supplementary Material ), water molecules further than 5.8Å from the nearest protein atom were also deleted, leaving 750 waters in a shell of approximately two- to three-layer thickness around the protein. The implicit hydration (GBSW) simulation is described in the Supplementary Material . In the production stage 25ns or 50ns trajectories were calculated. The CHARMM22/CMAP potential was used throughout.
Standard analyses of water dynamics and distributions were used to compare the PHS and GEO methods with the fully solvated protein simulation (PBC). Water motions are rapid, and so each trajectory frame was saved over a period of 100ps (2-fs intervals) after 25ns of simulation. Using the distance to the nearest protein surface atom, dnear, waters were classified as either belonging to the first shell (dnear<4Å), a second shell (4Å<dnear<7Å) or beyond reach. To capture the essence of water behavior in a certain shell, only water molecules were considered that had continuous residency in that particular shell (for 20ps in the first shell and 5ps in the second). Diffusion coefficients derived from mean-square displacements are given in Table 1. They are similar for the three simulations involving water, although the first-shell solvent molecules diffuse slightly faster in the PBC simulation. First-shell molecules have a smaller diffusion coefficient than those in the second shell, in agreement with previous studies 13,14.
The slower dynamics in the first shell is also obvious from Table 2, listing values from fitting the P2 rotational correlation function and the quantity
, the fraction of water molecules that still remain in a given shell after time t. The fitting function was
, except in the case of n(t) for the first shell, which was fitted to
. Here c accounts for the fraction of waters that stay in the first shell for a long time, presumably because they are either in the protein interior or bound to the surface. The number of these molecules was found to be 58 (
), 71 (
), and 55 (
) for PBC, PHS, and GEO, respectively. No fitting is reported here for the P1 correlation function as the decay timescale was longer than the analysis time intervals (20 and 5ps). The decay functions (shown in the Supplementary Material ), however, also suggest an acceptable agreement between the different models. Interestingly, ignoring the bound waters, exchange between the shells (mostly perpendicular to the protein surface) appears to be quite fast. This finding is, nevertheless, consistent with the diffusion and rotational correlation time as these are measured for waters that remain resident in a certain shell, preferentially sampling movement in parallel to the protein surface. In summary, we find that the solvent dynamics are overall very similar in the PHS and GEO calculations and close to those of the PBC simulation.
Distribution functions are a popular measure of solvent structure. It is known that the orientation of the water molecules close to a protein is not isotropic 13. In Fig. 1 we show the normalized distribution of two angles (θ and α); θ is the angle between the dipole moment of a water molecule and the vector connecting its oxygen to the center of geometry of the protein, indicating the orientation of waters with respect to the protein.
, and water-water orientation,
, for (a) first and (b) second shells.The distribution of water-protein angles, θ, shows small differences between the PBC and PHS/GEO simulations in the first shell. Second-shell waters show a slight orientational preference in the PBC simulation, favoring the oxygen to be closer to the protein (lysozyme has a net positive charge). A parallel or antiparallel dipole orientation, however, is not preferred for a noticeable fraction (∼10%) of second-shell waters in the PHS and GEO simulation. This behavior could arise from anisotropic hydrogen bonding that may exist at the solvent-vacuum boundary. Hybrid implicit/explicit solvation models have been developed to deal with this artifact (e.g., Lounnas et al. 8). We describe the orientation of a water molecule relative to others within 5Å by α, which is the angle between its dipole moment and those of neighboring waters belonging to the same shell. The distribution of α is similar in all cases and, as expected, is less anisotropic in the second shell as the influence of the protein surface decreases. Importantly, the orientation of waters toward each other does not appear to be affected by the solvent-vacuum boundary. Considering this and our analysis below, computationally costly approaches that randomize the dipole orientation of boundary waters may not be required.
For comparison of the structural and dynamic behavior of the protein several parameters were calculated over the last 10ns of the 25-ns trajectories. The main-chain root mean-square (RMS) deviations from the lysozyme crystal structure had equilibrated by then at 0.8, 1.4, 1.2, and 2.4Å for the PBC, PHS, GEO, and GBSW methods, respectively. Although the RMS deviations in the PHS and GEO calculations are somewhat greater than in the PBC simulation, longer (50-ns) simulations showed that these systems were stable. RMS fluctuations of Cα atoms in the shell and in the implicit simulation were compared to those in the PBC calculation. The resulting correlation coefficients (0.75, 0.80, and 0.71 for the PHS, GEO, and GBSW, respectively) showed a comparatively high agreement of main-chain dynamics.
Comparison with experimental data 11 is made in Table 3, listing correlation coefficients (R) and RMS differences (Δ) of the simulation-derived main-chain N-H and Asn/Gln side-chain NH2 Lipari-Szabo order parameters, S2 (data shown in Supplementary Material ), with the corresponding experimental values.
| Table 3 The correlation coefficients (R), and the RMS differences (Δ) between the simulation-derived and experimental order parameters |
| R (main) | Δ (main) | R (side) | Δ (side) | |||
|---|---|---|---|---|---|---|
| PBC | 0.67 | 0.060 | 0.74 | 0.155 | ||
| PHS | 0.62 | 0.061 | 0.75 | 0.160 | ||
| GEO | 0.53 | 0.074 | 0.72 | 0.157 | ||
| GBSW | 0.50 | 0.089 | 0.58 | 0.185 | ||
| See Buck et al. 11,12 for details on the derivation method. |
Although the main-chain correlation coefficients for the PBC and PHS are similar, the GBSW and GEO results are less correlated with the experimental data. The GEO method, however, shows a good agreement with the experiment for Asn/Gln side-chain dynamics. As noted before, the comparison between simulation and experimentally derived order parameters are not yet perfect 12. Remarkably the three methods that involve explicit waters give very similar results for the amplitude (S2) and also timescale (τe not shown) of the protein motions. This validates the inclusion of explicit solvent in simulations, but also suggests that principal solvation shell models reproduce the behavior seen in a fully solvated system.
In summary the data suggest that by comparison with full solvation under PBC, simple water shell models (PHS and GEO) result in acceptable water and protein behavior at a much lower computational cost. Using eight 3.2GHz Pentium 4 Xenon processors in parallel, 1.0ns of simulation took 36.5, 4.5, 3.4, and 11.3h for the PBC, PHS, GEO, and GBSW simulations, respectively. Water and protein behavior is almost identical between the water-shell and PBC approaches, whereas, as expected, there is a slight orientational preference of waters at the outer shell boundary. Overall, the results are encouraging for the further exploration of primary hydration-shell approaches in simulations of large proteins, protein docking, and protein-membrane systems.
We thank Mr. David Slochower for help with the early phases of the project.
The calculations were carried out on the High Performance Computing Cluster at Case Western Reserve University, supported with an award from the Provost’s office, and at the Ohio Supercomputing Center.
Note added in proof: The PHS and GEO simulations were extended to a longer period of 150ns, giving results closely similar to those reported.
1. (2002). Protein-water interactions in a dynamic world. Trends Biochem. Sci. 27, 203–208. Abstract | Full Text | PDF (388 kb) | CrossRef | PubMed
2. (2000). Solvent mobility and the protein ‘glass’ transition. Nat. Struct. Biol. 7, 34–38. CrossRef | PubMed
3. (1999). Effective energy function for proteins. Proteins 35, 133–152. CrossRef | PubMed
4. (2003). Generalized Born model with a simple smoothing function. J. Comput. Chem. 24, 1691–1702. CrossRef | PubMed
5. (1994). Dominant solvation effects from primary shell of hydration: approximation for molecular dynamics simulations. Biopolymers 35, 171–178. CrossRef | PubMed
6. (2000). Solvation in simulated annealing and high-temperature molecular dynamics of proteins: a restrained water droplet model. Int. J. Quantum Chem. 77, 174–186. PubMed
7. (1996). Hydrated myoglobin’s anharmonic fluctuations are not primarily due to dihedral transitions. Proc. Natl. Acad. Sci. USA 93, 55–59. CrossRef | PubMed
8. (1999). Towards molecular dynamics simulation of large proteins with a hydration shell at constant pressure. Biophys. Chem. 78, 157–182. CrossRef | PubMed
9. (2005). Evaluation of Poisson solvation models using a hybrid explicit/implicit solvent method. J. Phys. Chem. B 109, 5223–5236. PubMed
10. (2000). Application of the primary hydration shell approach to locally enhanced sampling simulated annealing: computer simulation of thyrotropin-releasing hormone in water. Biophys. J. 79, 66–79. Abstract | Full Text | PDF (1168 kb) | PubMed
11. (1995). Structural determinants of protein dynamics: analysis of 15N NMR relaxation measurements of main-chain and side-chain nuclei of hen egg white lysozyme. Biochemistry 34, 4041–4055. PubMed
12. (2006). Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme. Biophys. J. 90, L36–L38. Abstract | Full Text | PDF (111 kb) | CrossRef | PubMed
13. (2002). Molecular dynamics of water at the protein-solvent interface. J. Phys. Chem. B 106, 6617–6633. PubMed
14. (2006). Simulation studies of the protein-water interface. I. Properties at the molecular resolution. J. Chem. Phys. 124:234907, 1–18. PubMed