Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations

Abstract

Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 100-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of 2 and the all-hydrogen model. Continuum dielectric models of 1,1^*r, and 2^*r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreement than the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4^*r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.     

Protocol for MM/PBSA molecular dynamics simulations of proteins

Continuum solvent models have been employed in past years for understanding processes such as protein folding or biomolecular association. In the last decade, several attempts have been made to merge atomic detail molecular dynamics simulations with solvent continuum models. Among continuum models, the Poisson-Boltzmann solvent accessible surface area model is one of the oldest and most fundamental. Notwithstanding its wide usage for simulation of biomolecular electrostatic potential, the Poisson-Boltzmann equation has been very seldom used to obtain solvation forces for molecular dynamics simulation. We propose here a fast and reliable methodology to implement continuum forces in standard molecular mechanics and dynamics algorithms. Results for a totally unrestrained 1 ns molecular dynamics simulation of a small protein are quantitatively similar to results obtained by explicit solvent molecular dynamics simulations.     

Homology modelling and molecular dynamics simulations: comparative studies of human aquaporin-1

The structures of the mammalian water transport protein Aqp1 and of its bacterial homologue GlpF enables us to test whether homology models can be used to explore relationships between structure, dynamics and function in mammalian transport proteins. Molecular dynamics simulations (totalling almost 40 ns) were performed starting from: the X-ray structure of Aqp1; a homology model of Aqp1 based on the GlpF structure; and intermediate resolution structures of Aqp1 derived from electron microscopy. Comparisons of protein RMSDs vs. time suggest that the homology models are of comparable conformational stability to the X-ray structure, whereas the intermediate resolution structures exhibit significant conformation drift. For simulations based on the X-ray structure and on homology models, the flexibility profile vs. residue number correlates well with the crystallographic B-values for each residue. In the simulations based on intermediate resolution structures, mobility of the highly conserved NPA loops is substantially higher than in the simulations based on the X-ray structure or the homology models. Pore radius profiles remained relatively constant in the X-ray and homology model simulations but showed substantial fluctuations (reflecting the higher NPA loop mobility) in the intermediate resolution simulations. The orientation of the dipoles of water molecules within the pore is of key importance in maintaining low proton permeability through Aqp1. This property seems to be quite robust to the starting model used in the simulation. These simulations suggest that homology models based on bacterial homologues may be used to derive functionally relevant information on the structural dynamics of mammalian transport proteins.     

Implicit solvent models for flexible protein-protein docking by molecular dynamics simulation

The suitability of three implicit solvent models for flexible protein-protein docking by procedures using molecular dynamics simulation is investigated. The three models are (i) the generalized Born (GB) model implemented in the program AMBER6.0; (ii) a distance-dependent dielectric (DDD) model; and (iii) a surface area-dependent model that we have parameterized and call the NPSA model. This is a distance-dependent dielectric model modified by neutralizing the ionizable side-chains and adding a surface area-dependent solvation term. These solvent models were first tested in molecular dynamics simulations at 300 K of the native structures of barnase, barstar, segment B1 of protein G, and three WW domains. These protein structures display a range of secondary structure contents and stabilities. Then, to investigate the performance of the implicit solvent models in protein docking, molecular dynamics simulations of barnase/barstar complexation, as well as PIN1 WW domain/peptide complexation, were conducted, starting from separated unbound structures. The simulations show that the NPSA model has significant advantages over the DDD and GB models in maintaining the native structures of the proteins and providing more accurate docked complexes.     

Effective water model for Monte Carlo simulations of proteins

We present an effective theory for water. Our goal is to formulate on accurate model for the effects of solvation on protein dynamics, without incurring the huge computational cost and the slow temporal evolution typical of molecular dynamics simulations of liquids. We replace the individual water molecules in an all-atom potential with a local dielectric density field, with self interactions given by the Landau-Ginzburg free energy and external interactions by Lennard-Jones forces at the surface of the protein atoms. We explore conformational space with finite temperature Monte Carlo dynamics, using parallel Langevin and Fourier acceleration algorithms well suited to data-parallel computer architectures such as the Connection Machine. To establish the validity of our approximations, we compare our electrostatic contribution to the solvalion energy with the results of Lim, Bashford, and Karplus using a conventional static continuum dielectric cavity model, and the non electrostatic contributions with estimates of hydrophohic surface free energy. Our model can also accommodate ionic charges and temperature fluctuations, We propose future investigations extending our effective theory of solvation to include explicit orientational entropy and hydroxen-bonding terms. © 1995 John Wiley & Sons, Inc.     

H-bonding in protein hydration revisited

H-bonding between protein surface polar/charged groups and water is one of the key factors of protein hydration. Here, we introduce an Accessible Surface Area (ASA) model for computationally efficient estimation of a free energy of water-protein H-bonding at any given protein conformation. The free energy of water-protein H-bonds is estimated using empirical formulas describing probabilities of hydrogen bond formation that were derived from molecular dynamics simulations of water molecules at the surface of a small protein, Crambin, from the Abyssinian cabbage (Crambe abyssinica) seed. The results suggest that atomic solvation parameters (ASP) widely used in continuum hydration models might be dependent on ASA for polar/charged atoms under consideration. The predictions of the model are found to be in qualitative agreement with the available experimental data on model compounds. This model combines the computational speed of ASA potential, with the high resolution of more sophisticated solvation methods.     

Molecular dynamics simulations of aqueous pullulan oligomers

Molecular dynamics (MD) simulations have been used to model small-angle X-ray scattering (SAXS) data on aqueous solutions of four oligomeric segments of the glucan pullulan: the trimer G(3) (comprising one polymer repeating unit), the hexamer (G(3))(2), the nonamer (G(3))(3), and the dodecamer (G(3))(4). The AMBER force field was used in conjunction with the GB/SA continuum solvation model to calculate both the mean global dimensions of the oligomers from the limiting small angle scattering behavior and the shorter range structural information implicit in the Debye scattering function at larger scattering angles. This same force field and solvation treatment were employed earlier by Liu et al. (Macromolecules 1999, 32, 8611-8620) with apparent success in a rotational isomeric state (RIS) treatment of the same experimental data. The present work discloses that, despite numerical success in modeling the SAXS data, the RIS treatment, which includes only the interactions within dimeric segments of the polymer chain, fails to account accurately for excluded volume effects at the range of 3-12 sugar residues in the polymer backbone. It is suggested that MD simulations using continuum solvation models can be used to circumvent errors inherent in the computationally efficient RIS treatments of polymer nano- and picosecond dynamics while at the same time avoiding the heavy computational requirements of all-atom methods.     

Protein molecular dynamics with the generalized Born/ACE solvent model

Implicit solvent models are increasingly important for the study of proteins in aqueous solution. Here, the generalized Born (GB) solvent polarization model as implemented in the analytical ACE potential [Schaefer and Karplus (1996) J Phys Chem 100:1578] is used to perform molecular dynamics simulations of two small, homologous proteins: the immunoglobulin-binding domain of streptococcal protein G and the Ras binding domain of Raf. Several model parameterizations are compared through more than 60 ns of simulation. Results are compared with two simpler solvent models-an accessible surface area model and a distant-dependent dielectric model, with finite-difference Poisson calculations, with existing explicit solvent simulations, and with experimental data. The simpler models yield stable but distorted structures. The best GB/ACE implementation uses a set of atomic Voronoi volumes reported recently, obtained by averaging over a large database of crystallographic protein structures. A 20% reduction is applied to the volumes, compensating in an average sense for an excessive de-screening of individual charges inherent in the ACE self-energy and for an undersolvation of dipolar groups inherent in the GB screening function. This GB/ACE parameterization yields stable trajectories on the 0.5-1-ns time scale that deviate moderately (approximately 1.5-2.5 A) from the X-ray structure, reproduce approximately the surface distribution of charged, polar, and hydrophobic groups, and reproduce accurately backbone flexibility as measured by amide NMR-order parameters. Over longer time scales (1.5-3 ns), some of the protein G runs escape from the native energy basin and deviate strongly (3 A) from the native structure. The conformations sampled during the transition out of the native energy basin are overstabilized by the GB/ACE solvation model, as compared with a numerical treatment of the full dielectric continuum model.     

An Effective Solvation Term Based on Atomic Occupancies for Use in Protein Simulations

Abstract

Several approaches to the treatment of solvent effects based on continuum models are reviewed and a new method based on occupied atomic volumes (occupancies) is proposed and tested. The new method describes protein-water interactions in terms of atomic solvation parameters, which represent the solvation free energy per unit of volume. These parameters were determined for six different atoms types, using experimental free energies of solvation. The method was implemented in the GROMOS and PRESTO molecular simulation program suites. Simulations with the solvation term require 20-50% more CPU time than the corresponding vacuum simulations and are approximately 20 times faster than explicit water simulations. The method and parameters were tested by carrying out 200 ps simulations of BPTI in water, in vacuo, and with the solvation term. The performance of the solvation term was assessed by comparing the structures and energies from the solvation simulations with the equivalent quantities derived from several BPTI crystal structures and from the explicit water and vacuum simulations. The model structures were evaluated in terms of exposed total surface, buried and exposed polar surfaces, secondary structure preservation, number of hydrogen bonds, energy contributions, and positional deviations from BPTI crystal structures. Vacuum simulations produced unrealistic structures with respect to all criteria applied. The structures resulting from the simulations with explicit water were closer to the 5PTI crystal structure, although part of the secondary structure dissolved. The simulations with the effective solvation term produce structures that are normal according to all evaluations and in most respects are remarkably similar to the 5PTI crystal structure despite considerable positional fluctuations during the simulations. The segments where the model and crystal structures differ are known to be flexible and the observed difference may be physically realistic. The effective solvation term based on occupancies is not only very efficient in terms of computer time but also results in meaningful structural properties for BPTI. It may therefore be generally useful in molecular dynamics of macromolecules.     

Rapid boundary element solvation electrostatics calculations in folding simulations: successful folding of a 23-residue peptide.

Solvation effects play a profound role in the energetics of protein folding. While a continuum dielectric model of solvation may provide a sufficiently accurate estimate of the solvation effects, until now this model was too computationally expensive and unstable for folding simulations. Here we proposed a fast yet accurate and robust implementation of the boundary element solution of the Poisson equation, the REBEL algorithm. Using our earlier double-energy scheme, we included for the first time the mathematically rigorous continuous REBEL solvation term in our Biased Probability Monte Carlo (BPMC) simulations of the peptide folding. The free energy of a 23-residue beta beta alpha-peptide was then globally optimized with and without the solvation electrostatics contribution. An ensemble of beta beta alpha conformations was found at and near the global minimum of the energy function with the REBEL electrostatic solvation term. Much poorer correspondence to the native solution structure was found in the "control" simulations with a traditional method to account for solvation via a distance-dependent dielectric constant. Each simulation took less than 40 h (21 h without electrostatic solvation calculation) on a single Alpha 677 MHz CPU and involved more than 40,000 solvation energy evaluations. This work demonstrates for the first time that such a simulation can be performed in a realistic time frame. The proposed procedure may eliminate the energy evaluation accuracy bottleneck in folding simulations.     

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.     

Empirical solvation models in the context of conformational energy searches: application to bovine pancreatic trypsin inhibitor.

Continuum solvation models that estimate free energies of solvation as a function of solvent accessible surface area are computationally simple enough to be useful for predicting protein conformation. The behavior of three such solvation models has been examined by applying them to the minimization of the conformational energy of bovine pancreatic trypsin inhibitor. The models differ only with regard to how the constants of proportionality between free energy and surface area were derived. Each model was derived by fitting to experimentally measured equilibrium solution properties. For two models, the solution property was free energy of hydration. For the third, the property was NMR coupling constants. The purpose of this study is to determine the effect of applying these solvation models to the nonequilibrium conformations of a protein arising in the course of global searches for conformational energy minima. Two approaches were used: (1) local energy minimization of an ensemble of conformations similar to the equilibrium conformation and (2) global search trajectories using Monte Carlo plus minimization starting from a single conformation similar to the equilibrium conformation. For the two models derived from free energy measurements, it was found that both the global searches and local minimizations yielded conformations more similar to the X-ray crystallographic structures than did searches or local minimizations carried out in the absence of a solvation component of the conformational energy. The model derived from NMR coupling constants behaved similarly to the other models in the context of a global search trajectory. For one of the models derived from measured free energies of hydration, it was found that minimization of an ensemble of near-equilibrium conformations yielded a new ensemble in which the conformation most similar to the X-ray determined structure PTI4 had the lowest total free energy. Despite the simplicity of the continuum solvation models, the final conformation generated in the trajectories for each of the models exhibited some of the characteristics that have been reported for conformations obtained from molecular dynamics simulations in the presence of a bath of explicit water molecules. They have smaller root mean square (rms) deviations from the experimentally determined conformation, fewer incorrect hydrogen bonds, and slightly larger radii of gyration than do conformations derived from search trajectories carried out in the absence of solvent.     

Recent advances in implicit solvent-based methods for biomolecular simulations

Implicit solvent-based methods play an increasingly important role in molecular modeling of biomolecular structure and dynamics. Recent methodological developments have mainly focused on the extension of the generalized Born (GB) formalism for variable dielectric environments and accurate treatment of nonpolar solvation. Extensive efforts in parameterization of GB models and implicit solvent force fields have enabled ab initio simulation of protein folding to native or near-native structures. Another exciting area that has benefited from the advances in implicit solvent models is the development of constant pH molecular dynamics methods, which have recently been applied to the calculations of protein pK(a) values and the studies of pH-dependent peptide and protein folding.     

Unrestrained stochastic dynamics simulations of the UUCG tetraloop using an implicit solvation model.

Three unrestrained stochastic dynamics simulations have been carried out on the RNA hairpin GGAC[UUCG] GUCC, using the AMBER94 force field (Cornell et al., 1995. J. Am. Chem. Soc. 117:5179-5197) in MacroModel 5.5 (Mohamadi et al., 1990. J. Comp. Chem. 11:440-467) and either the GB/SA continuum solvation model (Still et al., 1990. J. Am. Chem. Soc. 112:6127-6129) or a linear distance-dependent dielectric (1/R) treatment. The linear distance-dependent treatment results in severe distortion of the nucleic acid structure, restriction of all hydroxyl dihedrals, and collapse of the counterion atmosphere over the course of a 5-ns simulation. An additional vacuum simulation without counterions shows somewhat improved behavior. In contrast, the two GB/SA simulations (1.149 and 3.060 ns in length) give average structures within 1.2 A of the initial NMR structure and in excellent agreement with results of an earlier explicit solvent simulation (Miller and Kollman, 1997. J. Mol. Biol. 270:436-450). In a 3-ns GB/SA simulation starting with the incorrect UUCG tetraloop structure (Cheong et al., 1990. Nature. 346:680-682), this loop conformation converts to the correct loop geometry (Allain and Varani, 1995. J. Mol. Biol. 250:333-353), suggesting enhanced sampling relative to the previous explicit solvent simulation. Thermodynamic effects of 2'-deoxyribose substitutions of loop nucleotides were experimentally determined and are found to correlate with the fraction of time the ribose 2'-OH is hydrogen bonded and the distribution of the hydroxyl dihedral is observed in the GB/SA simulations. The GB/SA simulations thus appear to faithfully represent structural features of the RNA without the computational expense of explicit solvent.     

Molecular dynamics studies of a DNA-binding protein: 2. An evaluation of implicit and explicit solvent models for the molecular dynamics simulation of the Escherichia coli trp repressor.

Although aqueous simulations with periodic boundary conditions more accurately describe protein dynamics than in vacuo simulations, these are computationally intensive for most proteins. Trp repressor dynamic simulations with a small water shell surrounding the starting model yield protein trajectories that are markedly improved over gas phase, yet computationally efficient. Explicit water in molecular dynamics simulations maintains surface exposure of protein hydrophilic atoms and burial of hydrophobic atoms by opposing the otherwise asymmetric protein-protein forces. This properly orients protein surface side chains, reduces protein fluctuations, and lowers the overall root mean square deviation from the crystal structure. For simulations with crystallographic waters only, a linear or sigmoidal distance-dependent dielectric yields a much better trajectory than does a constant dielectric model. As more water is added to the starting model, the differences between using distance-dependent and constant dielectric models becomes smaller, although the linear distance-dependent dielectric yields an average structure closer to the crystal structure than does a constant dielectric model. Multiplicative constants greater than one, for the linear distance-dependent dielectric simulations, produced trajectories that are progressively worse in describing trp repressor dynamics. Simulations of bovine pancreatic trypsin were used to ensure that the trp repressor results were not protein dependent and to explore the effect of the nonbonded cutoff on the distance-dependent and constant dielectric simulation models. The nonbonded cutoff markedly affected the constant but not distance-dependent dielectric bovine pancreatic trypsin inhibitor simulations. As with trp repressor, the distance-dependent dielectric model with a shell of water surrounding the protein produced a trajectory in better agreement with the crystal structure than a constant dielectric model, and the physical properties of the trajectory average structure, both with and without a nonbonded cutoff, were comparable.     

Ab initio folding of proteins with all-atom discrete molecular dynamics

Discrete molecular dynamics (DMD) is a rapid sampling method used in protein folding and aggregation studies. Until now, DMD was used to perform simulations of simplified protein models in conjunction with structure-based force fields. Here, we develop an all-atom protein model and a transferable force field featuring packing, solvation, and environment-dependent hydrogen bond interactions. We performed folding simulations of six small proteins (20-60 residues) with distinct native structures by the replica exchange method. In all cases, native or near-native states were reached in simulations. For three small proteins, multiple folding transitions are observed, and the computationally characterized thermodynamics are in qualitative agreement with experiments. The predictive power of all-atom DMD highlights the importance of environment-dependent hydrogen bond interactions in modeling protein folding. The developed approach can be used for accurate and rapid sampling of conformational spaces of proteins and protein-protein complexes and applied to protein engineering and design of protein-protein interactions.     

Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

A continuum electrostatics approach for molecular dynamics (MD) simulations of macromolecules is presented and analyzed for its performance on a peptide and a globular protein. The approach incorporates the screened Coulomb potential (SCP) continuum model of electrostatics, which was reported earlier. The model was validated in a broad set of tests some of which were based on Monte Carlo simulations that included single amino acids, peptides, and proteins. The implementation for large-scale MD simulations presented in this article is based on a pairwise potential that makes the electrostatic model suitable for fast analytical calculation of forces. To assess the suitability of the approach, a preliminary validation is conducted, which consists of (i) a 3-ns MD simulation of the immunoglobulin-binding domain of streptococcal protein G, a 56-residue globular protein and (ii) a 3-ns simulation of Dynorphin, a biological peptide of 17 amino acids. In both cases, the results are compared with those obtained from MD simulations using explicit water (EW) molecules in an all-atom representation. The initial structure of Dynorphin was assumed to be an alpha-helix between residues 1 and 9 as suggested from NMR measurements in micelles. The results obtained in the MD simulations show that the helical structure collapses early in the simulation, a behavior observed in the EW simulation and consistent with spectroscopic data that suggest that the peptide may adopt mainly an extended conformation in water. The dynamics of protein G calculated with the SCP implicit solvent model (SCP-ISM) reveals a stable structure that conserves all the elements of secondary structure throughout the entire simulation time. The average structures calculated from the trajectories with the implicit and explicit solvent models had a cRMSD of 1.1 A, whereas each average structure had a cRMSD of about 0.8A with respect to the X-ray structure. The main conformational differences of the average structures with respect to the crystal structure occur in the loop involving residues 8-14. Despite the overall similarity of the simulated dynamics with EW and SCP models, fluctuations of side-chains are larger when the implicit solvent is used, especially in solvent exposed side-chains. The MD simulation of Dynorphin was extended to 40 ns to study its behavior in an aqueous environment. This long simulation showed that the peptide has a tendency to form an alpha-helical structure in water, but the stabilization free energy is too weak, resulting in frequent interconversions between random and helical conformations during the simulation time. The results reported here suggest that the SCP implicit solvent model is adequate to describe electrostatic effects in MD simulation of both peptides and proteins using the same set of parameters. It is suggested that the present approach could form the basis for the development of a reliable and general continuum approach for use in molecular biology, and directions are outlined for attaining this long-term goal.     

Temperature-dependent molecular dynamics and restrained X-ray refinement simulations of a Z-DNA hexamer

Molecular dynamics simulations of the Z-DNA hexamer 5BrdC-dG-5BrdC-dG-5BrdC-dG were performed at several temperatures between 100 K and 300 K. Above 250 K, a strong sequence-dependent flexibility in the nucleic acid is observed, with the guanine sugar and the phosphate of GpC sequences much more mobile than the cytosine sugar and phosphate of CpG sequences. At 300 K, the hexamer is in dynamic equilibrium between several Z forms, including the crystallographically determined ZI and ZII forms. The local base-pair geometry, however, is not very variable, except for the roll of the base-pairs. The hexamer molecular dynamics trajectories have been used to test the restrained parameter crystallographic refinement model for nucleic acids. X-ray diffraction intensities corresponding to observed diffraction data were computed. The average structures obtained from the simulations were then refined against the calculated intensities, using a restrained least-squares program developed for nucleic acids in order to analyse the effects of the refinement model on the derived quantities. In general, the temperature dependence of the atomic fluctuations determined directly from the refined Debye-Waller factors is in reasonably good agreement with the results obtained by calculating the atomic fluctuations directly from the Z-DNA molecular dynamics trajectories. The agreement is best for refinement of temperature factors without restraints. At the highest temperature studied (300 K), the effect of the refinement on the most mobile atoms (phosphates) is to significantly reduce the mean-square atomic fluctuations estimated from the refined Debye-Waller factors below the actual values (less than (delta r)2 greater than congruent to 0.5 A2). Analysis of the temperature-dependence of the mean-square atomic fluctuations provides information concerning the conformational potential within which the atoms move. The calculated temperature-dependence and anharmonicity of the Z-DNA helix are compared with the results observed for proteins. The average structures from the simulations were refined against the experimental X-ray intensities. It is found that low-temperature molecular dynamics simulations provide a useful tool for optimizing the refinement of X-ray structures.     

Acceptable protein and solvent behavior in primary hydration shell simulations of hen lysozyme

The "primary hydration shell" method in molecular dynamics simulations uses a two- to three-layer thick shell of explicitly represented water molecules as the solvent around the protein of interest. We show that despite its simplicity, this computationally cheap model is capable of predicting acceptable water and protein behavior using the CHARMM22/CMAP potential function. For protein dynamics, comparisons are made with Lipari-Szabo order parameters. These have been derived from NMR relaxation parameters for pico-nano second motions of the NH groups in the main-chain and NH(2) groups in Asn/Gln side chains in hen lysozyme. It is also shown that an even simpler, and therefore faster, water-shell model leads to results in similarly good agreement with experiments, and also compared with simulations using a full box of water with periodic boundary conditions or with an implicit solvation model. Thus, the primary hydration shell method should be useful in making larger systems accessible to extensive simulations.