Force Field Benchmark of Amino Acids: I ... - ACS Publications

Article Cite This: J. Chem. Inf. Model. 2018, 58, 1037−1052

pubs.acs.org/jcim

Force Field Benchmark of Amino Acids: I. Hydration and Diffusion in Different Water Models Haiyang Zhang,† Chunhua Yin,† Yang Jiang,‡ and David van der Spoel*,§

Downloaded via KAOHSIUNG MEDICAL UNIV on November 15, 2018 at 21:50:38 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

Department of Biological Science and Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, Beijing 100083, China ‡ Beijing Key Lab of Bioprocess, College of Life Science and Technology, Beijing University of Chemical Technology, Box 53, Beijing 100029, China § Uppsala Center for Computational Chemistry, Science for Life Laboratory, Department of Cell and Molecular Biology, Uppsala University, Husargatan 3, Box 596, SE-75124 Uppsala, Sweden S Supporting Information *

ABSTRACT: Thermodynamic and kinetic properties are of critical importance for the applicability of computational models to biomolecules such as proteins. Here we present an extensive evaluation of the Amber ff99SB-ILDN force field for modeling of hydration and diffusion of amino acids with three-site (SPC, SPC/E, SPC/Eb, and TIP3P), four-site (TIP4P, TIP4P-Ew, and TIP4P/2005), and five-site (TIP5P and TIP5P-Ew) water models. Hydration free energies (HFEs) of neutral amino acid side chain analogues have little dependence on the water model, with a root-mean-square error (RMSE) of ∼1 kcal/mol from experimental observations. On the basis of the number of interacting sites in the water model, HFEs of charged side chains can be putatively classified into three groups, of which the group of three-site models lies between those of four- and five-site water models; for each group, the water model dependence is greatly eliminated when the solvent Galvani potential is considered. Some discrepancies in the location of the first hydration peak (RRDF) in the ion−water radial distribution function between experimental and calculated observations were detected, such as a systematic underestimation of the acetate (Asp side chain) ion. The RMSE of calculated diffusion coefficients of amino acids from experiment increases linearly with the increasing diffusion coefficients of the solvent water models at infinite dilution. TIP3P has the fastest diffusivity, in line with literature findings, while the “FB” and “OPC” water model families as well as TIP4P/2005 perform well, within a relative error of 5%, and TIP4P/2005 yields the most accurate estimate for the water diffusion coefficient. All of the tested water models overestimate amino acid diffusion coefficients by approximately 40% (TIP4P/2005) to 200% (TIP3P). Scaling of protein−water interactions with TIP4P/ 2005 in the Amber ff99SBws and ff03ws force fields leads to more negative HFEs but has little influence on the diffusion of amino acids. The most recent FF/water combinations of ff14SB/OPC3, ff15ipq/SPC/Eb, and fb15/TIP3P-FB do not show obvious improvements in accuracy for the tested quantities. These findings here establish a benchmark that may aid in the development and improvement of classical force fields to accurately model protein dynamics and thermodynamics. properties of water molecules as accurately as possible.3−8 A special Web site with information on the science of water is maintained by Martin Chaplin.9 To date, however, many of the

1. INTRODUCTION Extensive attention has been paid to the study of molecular interactions in water using both experimental and computational tools, attempting to explore the role of water in biology.1,2 In the past decades, a large number of theoretical potentials have been proposed to capture the physical © 2018 American Chemical Society

Received: January 19, 2018 Published: April 12, 2018 1037

DOI: 10.1021/acs.jcim.8b00026 J. Chem. Inf. Model. 2018, 58, 1037−1052

Article

Journal of Chemical Information and Modeling

advantageous for modeling properties of disordered proteins.37−40 The latest Amber protein FFs, ff15ipq18 and fb15,41 have been built for consistency with the recent highquality water models SPC/Eb19 and TIP3P-FB,28 and shown improved performance for diffusion and temperature-dependent properties of proteins, respectively. Despite the quality of water models, a key observation in recent years has been the imbalance between protein−water and protein−protein interactions. Simulated disordered peptides were found to be structurally too compact relative to experiment,20,42,43 and artificial protein aggregation was observed under self-crowding conditions in simulations.44−46 These observations indicate that protein−water interactions are likely too weak compared with protein−protein interactions, as confirmed by a slightly more positive hydration free energy of amino acid side chain analogues.20,47,48 To restore the balance, some researchers have taken a first step and proposed a simple strategy to tune short-range solute−water Lennard-Jones (LJ) interactions, allowing for predictions that are more consistent with experimental solvation free energies49 or small-angle X-ray scattering and Förster resonance energy transfer measurements.20 To capture thermodynamic and kinetic properties of biomolecules like proteins in water accurately, a high-quality water model and balanced interactions between interacting molecules are therefore required, and more importantly, thorough identification of the origin of force field deficiencies such as hydration and diffusion properties is highly desirable. Here we present an extensive benchmark of hydration and diffusivity of amino acid analogues for the Amber protein force field ff99SB-ILDN10 paired with nine water models, namely, SPC,14 SPC/E,15 SPC/Eb,19 TIP3P,13 TIP4P,13 TIP4P-Ew,25 TIP4P/2005,27 TIP5P,16 and TIP5P-Ew.26 For TIP4P/2005, additional variants of Amber FFs, including ff99SBw,20 ff99SBws,20 ff03w,33 and ff03ws,20 were examined as well. The latest Amber FFs, ff14SB,50 ff15ipq,18 and fb15,41 were tested with the recent water models OPC3,30 SPC/Eb,19 and TIP3P-FB,28 respectively. The similarities, differences, and/or deficiencies of force field/water combinations are discussed in detail. This work can be useful for further efforts in protein force field development.

popular models reproduce water properties with limited accuracy only, and thus, the choice of water model in molecular dynamics (MD) simulations commonly hinges on the research purpose and the compatibility with the used force fields. For instance, the native water models of the Amber,10 GROMOS,11 and OPLS12 force fields (FFs) are TIP3P (transferable intermolecular potential 3P),13 SPC (simple point charge),14 and TIP4P,13 respectively; more water models are supported as well for Amber, such as SPC/E (extended SPC),15 TIP4P,13 and TIP5P.16 Of these, SPC/E appeared to be compatible with varied biomolecular FFs.17 It could be argued, however, that SPC/E may be incompatible with protein force fields because the polarization energy is built into its parametrization.15 No biomolecular force fields were designed with an equivalent approach, although the recent Amber force field ff15ipq18 was tuned with reference to the SPC/Eb19 cousin water model. Most of the classical additive force fields like Amber largely followed historical precedent and were developed and tested with an earlier generation of water models, such as the TIP3P model developed 35 years ago.13 Clearly, one of the weaknesses of these FFs is the reliance on the water model used.20 A wellknown drawback of TIP3P is its fast self-diffusivity and low viscosity,21−23 indicating that TIP3P is not suitable for predicting absolute kinetic properties of biomolecules in water. For this reason, diffusion/viscosity corrections are required for a meaningful comparison with experimental kinetic properties in studies of, for instance, metal ions and protein folding dynamics.23,24 Refined water models were proposed a decade ago that were carefully parametrized against extensive experimental data, for instance, TIP4P-Ew (adaption of TIP4P for use with Ewald techniques),25 TIP5P-Ew (adaption of TIP5P for Ewald),26 and TIP4P/2005.27 The SPC/Eb model, with a slight increase of 1% in the O−H bond length compared with the original SPC/E model, was proposed recently in order to obtain improved reproduction of rotational diffusion and NMR spectral density of pure water and proteins.19 Two water model families of the latest generation, “FB” (like TIP3P-FB and TIP4P-FB)28 and “OPC” (like OPC3 and OPC),29,30 were developed using different strategies such as force balance and global optimization, respectively. These models were shown to reproduce a comprehensive set of liquid bulk properties accurately. However, it is difficult to choose a water model for biomolecular simulations without considering the biomolecular force field.31 Using an improved water model should, in principle, result in enhanced accuracy of the simulated properties of biomolecules, since the hydrophobic effect that is the driving force for biomolecular associations and folding is largely determined by solvent properties.32 By introducing modifications for protein FFs, several groups have made progress for application of these FFs with the recently optimized water models.20,33,34 With minimal backbone changes, Best and Mittal33 combined Amber ff03 with TIP4P/2005 (the resulting FF was denoted as ff03w) to reproduce helix propensity; Nerenberg and Head-Gordon34 adapted Amber ff99SB35 to TIP4P-Ew in order to obtain improved conformational ensembles of small peptides. Both of these adapted FFs show better performance than the native water model TIP3P in combination with the unchanged Amber FF. Similarly, for a faithful reproduction of helix propensity by ff99SB-ILDN,10 Best and co-workers then adjusted backbone torsions and refitted atomic charges of charged amino acids and side chain torsions of Asp for use with TIP4P/2005 (ff99SBw).20,36 These adaptions have been shown to be

2. COMPUTATIONAL METHODS 2.1. Structural Model. All of the amino acid residues (except Gly and Pro) that are available in the Amber force fields were examined in this work. Construction of side chain analogues followed Shirts and co-workers,48 and the net charges for neutral and charged analogues were forced to be zero and unity (−1 for Asp and Glu and +1 for Arg, His, and Lys), respectively, by adjustment of the atomic charge of the βcarbon (here called the method of Cβ modification).51 Two tautomeric states of neutral His, namely, Hid and Hie, with a proton connected to the δ- or ε-nitrogen, respectively, were constructed separately. Similarly, three partial analogues of guanidinium (Gud+), imidazolium (Imd+), and ammonium (NH4+) ions were also built on the basis of Arg+, His+, and Lys+, respectively, for direct comparison with experimental observations. These side chain analogues were used to calculate hydration properties of amino acids. Experimental diffusion constants of amino acids at infinite dilution used here for reference were determined in a buffer of pH ∼3.5,52 at which amino acids usually exist in zwitterionic form (NH3+−X−COO−). The original Amber FFs do not have parameters for zwitterionic amino acid residues, and we 1038


Article


0.002 ps. Electrostatic interactions were handled by the particle mesh Ewald (PME) method60,61 with a real-space cutoff of 1.0 nm, and Lennard−Jones interactions were truncated at 1.0 nm. The velocity rescaling thermostat62 and Parrinello−Rahman barostat63,64 were applied to maintain the temperature and pressure at 298.15 K and 1 bar with coupling time constants of 1 and 5 ps, respectively. The SETTLE algorithm65 was used to maintain a rigid water model. The calculation of diffusion coefficients should ideally be done on the basis of simulations in the NVE ensemble. Because of the practical issues of NVE simulations (see more discussion in the work by York and coworkers,23 who concluded that the computed diffusion coefficients of water and Mg2+ and the temperature distribution in NVE simulations are very close to those in NVT simulations), we chose the NVT ensemble for the calculations. After an energy minimization and 1 ns NVT, 1 ns NPT, and 3 ns NVT equilibrations, production trajectories were performed at NVT for 20 ns with atomic coordinates saved to disk every 2 ps. Finite-size effects under periodic boundary conditions (PBCs) have been shown to influence the determination of the translational diffusion constant from simulations for molecules such as water and metal ions.21,23,66,67 In order to correct for such effects, the box-size-dependent diffusion constants (DPBC) were used to extrapolate the diffusion constant to the infinite dilution limit (D0), as shown in eq 2:21,68

constructed them manually on the basis of the available N- and C-terminal amino acids. The carbonyl group of N-terminal amino acid building blocks was changed to a carboxylate group (−COO−), of which atom types were assigned in accordance with the charged Asp residue, resulting in zwitterionic amino acids. Atomic charges for carboxylate oxygen atoms were set to the average value for all of the carboxylate groups of C-terminal amino acids. Finally, if necessary, the charge of the carbon atom of −COO− was adjusted to ensure an integral net charge. This manual assignment is called the method of CO modification. These zwitterionic forms were used to evaluate the diffusion properties of amino acids. 2.2. Amino Acid Force Field and Water Models. The Amber ff99SB-ILDN force field10 was used to model side chain analogues and zwitterionic forms of amino acids in combination with three catalogues of molecular mechanics water models, namely, three-site (SPC,14 SPC/E,15 SPC/Eb,19 and TIP3P13), four-site (TIP4P,13 TIP4P-Ew,25 and TIP4P/200527), and fivesite (TIP5P16 and TIP5P-Ew26) models. Four other variants of the Amber protein FF, namely, ff99SBw,20 ff99SBws,20 ff03w,33 and ff03ws,20 were examined as well. These FFs were adapted for use with the optimized TIP4P/2005 water model27 by Best and co-workers, and the resulting FFs showed good performances for the properties of helix propensity, disordered proteins, and nonspecific protein associations in protein simulations.20,33,36,53 Both ff99SBw and ff99SBws were based on ff99SB*-ILDN-Q,36 while the latter used a scaling of protein−water LJ interactions (indicated by the suffix “s”), as computed by

εOi = αεOLBi = α εOεi

DPBC = D0 −

2.83729kBT 6πηL

(2)

(1)

where kB is Boltzmann’s constant, T is the absolute temperature (298.15 K), η is the shear viscosity of the solvent, and L is the box length. In practice, the y intercept of a linear fit of a plot of DPBC versus 1/L gives D0, and the corresponding slope can be used to calculate η. It should be noted that the viscosity was shown to be system-size independent; given a solvent viscosity, eq 2 can also be used to correct the system-size dependences of diffusion coefficient calculations under PBCs.21,69 Similar to the method used by York and co-workers,23 each 20 ns NVT production trajectory was divided into 100 blocks equally. For each block, mean squared displacements (MSDs) of the diffusing particle were calculated and averaged over all of the molecules, and the box-size-dependent diffusion constant was then calculated from the MSD using the Einstein relation.70 The final DPBC for each box was obtained by averaging over all of the blocks, and error estimates of DPBC were evaluated by block averaging.69 With the aid of the R program,71 D0 and η and their corresponding errors were determined from linear regression of DPBC versus 1/L (eq 2) by generating 1000 data sets from a normal distribution of DPBC randomly and obtaining the corresponding averages and standard deviations. 2.3.2. Dielectric Permittivity and Multipole Moment. Static relative dielectric constants (ε0) at 298.15 K and 1 bar were computed from the fluctuations of the total dipole moment of the simulation box.72,73 The molecular dipole moment (μ) and dielectric constant (ε0) of water were calculated with the GROMACS analysis tool “gmx dipoles”.56−59 The quadrupole moments (QT) of three- and four-site water models were calculated following the work by Abascal and Vega (eq 5 in ref 74), and that of the five-site models was obtained from ref 16. 2.3.3. Surface Potential and Surface Tension. After the 1 ns NPT equilibration, a 3.1 nm × 3.1 nm × 9.3 nm slab with 1000

where εO and εi are LJ parameters (ε) for water oxygen and protein atom i, εOi is the result for mixed atom types generated by the Lorentz−Berthelot (LB) combination rule, and α is the scaling factor (α = 1.1 for ff99SBws and ff03ws).20 Similarly, both ff03w and ff03ws were based on ff0354 with a slight backbone modification,53 and the latter implemented scaled protein−water interactions similar to ff99SBws. It should be noted that the atomic charges of neutral, N-, and C-terminal amino acids remained unchanged in ff99SBw(s)20 and were identical to those in ff99SB-ILDN,10 while ff03w(s)20,33 shared identical atomic charges with ff0354 for all amino acid residues. FF/water combinations of latest Amber FFs and water models, ff14SB/OPC3, ff15ipq/SPC/Eb, and fb15/TIP3P-FB (extracted from AmberTools 1655), were used to model side chain analogues and zwitterionic amino acid residues as well. It should be noted that the Amber ff14SB is not built around OPC3, and because of the too-fast self-diffusivity of TIP3P,21−23 we choose OPC3 to examine whether the diffusion properties of amino acids can be improved by using the latest water models. Water properties for two more four-site models, TIP4P-FB28 and OPC,29 were also calculated for a comparison with the other models. 2.3. Estimation of Water Properties. 2.3.1. Diffusion Constants. All of the MD simulations were carried out with full periodic boundary conditions using the GROMACS package (version 5.1).56−59 For estimation of pure-water properties, seven cubic boxes containing 512, 768, 1024, 1536, 2048, 3072, and 4096 pure water molecules with box lengths of about 2.5, 2.9, 3.2, 3.6, 4.0, 4.6, and 5.0 nm, respectively, were simulated separately. The leapfrog algorithm was used to integrate classical Newton’s equations of motion using a time step of 1039


Article

Journal of Chemical Information and Modeling pure water molecules was simulated at NVT for 10 ns with a protocol close to that in section 2.3.1. The last 8 ns was used for the calculation of the surface potential (χ) and surface tension (γ). The GROMACS tool “gmx potential” was used to compute the potential across the box, where the periodic boundaries were taken into account to make the potential periodic as well. The potential was computed by integrating twice the charge density per slab, and it was set to zero in the gas phase. A script was used to compute the average potential in the water phase. By means of the same simulations, the surface tension for the water models was computed from the difference in the lateral and normal pressures (see ref 75 for details). 2.4. Estimation of Amino Acid Hydration. Hydration free energies (HFEs) of amino acid side chain analogues were calculated by thermodynamic integration (TI). A Langevin (stochastic dynamics) integrator76 was used for TI simulations, with a friction coefficient of 1 ps−1. The LINCS algorithm77 was used to constrain all bonds of the solute molecules. A two-stage TI was performed along an alchemical reaction coordinate of λ = 0, 0.25, 0.5, 0.75, and 1 for decoupling of Coulombic interactions and then λ = 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, and 1 for decoupling of van der Waals interactions. Each λ was simulated for 2 ns at 298.15 K and 1 bar, and the first 100 ps was removed for system equilibration. Each system contains only one side chain analogue and 640 water molecules in a cubic box. The protocol for TI simulations was similar to that for the above pure-water calculations, and more details can be found in ref 78. HFE calculations that are typically done under PBCs (as in this work) do not have any boundary between air and the liquid and are known as intrinsic free energies (ΔGintr). A surface contribution (ΔGsurf) from crossing the vacuum−bulk interface must be added, giving rise to the physical (real) free energy of solvation (eq 3):79−81 ΔGreal = ΔGintr + ΔGsurf = ΔGintr + qØG

molecules. To monitor water structures around amino acids for clarity, Here an unnormalized RDF was used: ρ(r) = n(r)/ (4πr2Δr), where r is the distance between the reference site of an ion and the oxygen atom of a water molecule and n(r) is the average number of the given particle in a spherical shell with thickness Δr at a distance r.86 ρ(r) is also known as the water number density (nm−3) around ions at r. The location of the first hydration peak (RRDF) in ρ(r) can be used for comparison with experimental estimates if available.87−91 Nitrogen and oxygen atoms were chosen as the reference sites for cations and anions, respectively.11 2.5. Estimation of Amino Acid Diffusion. Similar to the water diffusion calculations, each zwitterionic amino acid was immersed in a cubic box containing 512, 768, 1024, 2048, or 4096 water molecules separately. After energy minimization and equilibrations, 30 separate simulations were performed for each box and each amino acid, and these were used to determine the box-size-dependent diffusion constants (DPBC). Linear fits of plots of DPBC versus 1/L (eq 2) give the calculated translational diffusion coefficients (Dcalc 0,AA) of the amino acids at infinite dilution by extrapolation. In order to correct for potential errors from the diffusion of water models, a scaled 23 quantity (Dscaled 0,AA ) was obtained using eq 4: scaled calc D0,AA = D0,AA

exp D0,W calc D0,W

(4)

calc where Dexp 0,W and D0,W are experimental and calculated water diffusion coefficients, respectively. This correction allows a more meaningful comparison with experimental observations.

3. RESULTS 3.1. Water Properties. Simulated self-diffusion coefficients (DPBC) of water molecules as functions of the inverse box length (1/L) examined are presented in Figure 1. As expected, the water diffusion constants depend on the simulated system size linearly and increase as the box length increases. The extrapolated diffusion constants at infinite dilution (D0) and the corresponding viscosities computed using eq 2 are summarized in Table 1. Among the 13 water models, TIP3P gives the largest diffusion coefficient, D0 = (6.14 ± 0.05) × 10−5 cm2/s, which significantly overestimates the experimental determination (D0 = 2.299 × 10−5 cm2/s),92 and TIP4P/2005 predicts a value of D0 = (2.31 ± 0.02) × 10−5 cm2/s, which reproduces the experiment perfectly (Figure 1a). The latest two families of water models, “FB” and “OPC”, yield predictions close to the experimental value, yet with a slight underestimation and overestimation by about 5%, respectively (Figure 1b and Table 1). SPC/Eb gives the smallest diffusion coefficient at infinite dilution, and the D0 values (in units of 10−5 cm2/s) of the water models decrease in the order TIP3P (6.14) > SPC (4.52) > TIP4P (3.86) > TIP5P-Ew (3.09) > TIP5P (2.93) > SPC/E (2.76) > TIP4P-Ew (2.63) > OPC3 (2.41) > OPC (2.37) > TIP4P/2005 (2.31) > TIP4P-FB (2.21) > TIP3P-FB (2.20) > SPC/Eb (2.10) (Figure 1 and Table 1). Our D0 calculations are in excellent agreement with the results published by others (with the finite-size effects under PBCs corrected) for TIP3P (6.05 × 10−5 cm2/s) and SPC/E (2.79 × 10−5 cm2/s)21,23 and are slightly smaller than the predictions for TIP4P-Ew (2.72 × 10−5 cm2/s),23 SPC/E (2.97 × 10−5 cm2/s), and TIP4P/2005 (2.49 × 10−5 cm2/s).67 The original SPC/Eb work reported a value of 2.34 × 10−5 cm2/s in the NPT ensemble,19, which is 11% larger than our prediction in the NVT ensemble.

(3)

where q is the net charge of the solute and ØG is the solvent Galvani potential jump in the air to liquid direction (ØG = Fχ, where F = is the Faraday constant and χ is the interface potential). For charged side chains, finite-size artifacts and potential inaccuracies in the solvent model permittivity need to be considered with care as well.82 It was noted that these two corrections are quite small and negligible for PME/PBC simulations (as we did in this work).78,83 Experimental observations for the HFEs of ions are based on a debated hydration free energy of the proton, ΔGhyd * (H+), ranging from −265 to −250 kcal/mol,84 and to date there is no consensus on the true value. Our previous studies on charged amino acids and molecular ions78,83 indicated that the calculated real free energy of ionic hydration (i.e, taking the solvent Galvani potential into account) with classical molecular mechanics force fields appeared to best match the experiments using a ΔGhyd * (H+) of −262.4 kcal/mol85 for reference. This value was therefore used here for comparison with the calculations. Besides the solvation free energies, the water structure around amino acids is of crucial importance for protein dynamics in water. Charged side chain analogues of Arg+, His+, Imd+, Lys+, NH4+, Asp−, and Glu− were simulated for 10 ns separately without counterions in a cubic box of 640 water molecules, and the last 8 ns was used to compute radial distribution functions (RDFs) of a given ion around water 1040


Article


water models show more or less acceptable performance for the water viscosity (Table 1). It should be noted that reverse nonequilibrium molecular dynamics (RNEMD), a method close to the real experimental technique, produced dramatically small estimates of η = 0.33 mPa·s for SPC/E and 0.58 mPa·s for TIP4P/2005, whereas for SPC, TIP3P, TIP4P, TIP4P-Ew, TIP5P, and TIP5P-Ew comparable results were obtained by RNEMD.99 The viscosities of the water models examined are found to correlate strongly (R = −0.96) with their self-diffusion coefficients (D0), and η deceases linearly with increasing D0, as indicated by the blue circles in Figure 2.

Figure 1. Box-size-dependent diffusion constants of water calculated from MD simulations with varied box lengths (L) under periodic boundary conditions (PBCs) for (a) early-generation and (b) latest water models. Values for TIP4P/2005 are given in both panels for comparison. The solid lines indicate linear fits of data points; the intercept (i.e., where L approaches infinity) gives the diffusion constant at infinite dilution, and the slope is used to estimate the viscosity of the water model (eq 2).

Figure 2. Correlations of calculated water diffusion coefficients with the corresponding shear viscosities of the solvent water (blue circles) and with the root-mean-square errors (RMSEs) (red triangles) of the calculated amino acid diffusion coefficients from experimental observations at infinite dilution. Solid lines indicate linear fits of data points. R is the correlation coefficient, and RMSE is for the calculations with ff99SB-ILDN.

Linear fits of DPBC vs 1/L (Figure 1 and eq 2) produce estimates for the shear viscosities (η) of the water models, which are listed in Table 1. Compared with the experimental value (η = 0.89 mPa·s),95,96 the TIP3P, SPC, and TIP4P models severely underestimate the water viscosity, and among these, TIP3P yields the worst (smallest) viscosity of 0.32 ± 0.04 mPa·s, in line with reports by others.21,98−102 TIP4P/2005 reproduces the viscosity well (η = 0.84 ± 0.08 mPa·s), as confirmed by several groups.67,98,100,101 OPC3, SPC/Eb, and TIP4P-FB show good performance for η as well, and the other

Surface tensions are underestimated to some extent for all of the water models, of which the four-site OPC model performs best and predicts a surface tension of 68.2 mN/m, which is very close to the experimental value of 71.99 mN/m.97 This underestimation can be ascribed to the lack of polarization in the models and omission of long-range dispersion.103,104 Compared with the experimental dipole moment of liquid

Table 1. Calculated Properties of the Water Models Tested at 298.15 Ka OPC3 SPC SPC/E SPC/Eb TIP3P TIP3P-FB OPC TIP4P TIP4P-Ew TIP4P-FB TIP4P/2005 TIP5P TIP5P-Ew expt

μ (D)

QT (10−1 D nm)

ε0

D0,W (10−5 cm2/s)

η (mPa·s)

χ (V)

ØG (kcal/mol/e)

γ (mN/m)

2.430 2.274 2.351 2.374 2.347 2.419 2.480 2.178 2.321 2.429 2.305 2.292 2.292 2.9 (0.6)b

2.060 1.969 2.036 2.076 1.721 2.052 2.300 2.147 2.164 2.151 2.297 1.565 1.565

77.8 (1.9) 65.1 (1.3) 70.6 (1.6) 70.5 (1.3) 96.2 (1.8) 78.6 (1.8) 76.3 (1.9) 51.1 (0.8) 63.4 (1.4) 76.6 (2.1) 56.7 (1.6) 89.9 (3.7) 94.4 (3.6) 78.375c

2.41 (0.02) 4.52 (0.03) 2.76 (0.02) 2.10 (0.01) 6.14 (0.05) 2.20 (0.01) 2.37 (0.02) 3.86 (0.03) 2.63 (0.04) 2.21 (0.01) 2.31 (0.02) 2.93 (0.03) 3.09 (0.03) 2.299d

0.87 (0.07) 0.45 (0.03) 0.78 (0.07) 0.91 (0.06) 0.32 (0.04) 1.00 (0.08) 1.01 (0.09) 0.55 (0.07) 0.75 (0.07) 0.90 (0.06) 0.84 (0.08) 0.82 (0.09) 0.66 (0.07) 0.89e

−0.624 −0.605 −0.589 −0.590 −0.534 −0.569 −0.544 −0.523 −0.513 −0.517 −0.486 −0.104 −0.102

−14.4 −14.0 −13.6 −13.6 −12.3 −13.1 −12.5 −12.1 −11.8 −11.9 −11.2 −2.4 −2.4

60.0 (0.2) 49.5 (0.3) 56.8 (0.4) 60.4 (0.3) 45.9 (0.7) 59.8 (0.3) 68.2 (0.2) 51.4 (0.6) 58.0 (0.2) 63.7 (0.3) 62.6 (0.6) 47.9 (0.6) 51.5 (0.4) 71.99f

a

Definitions: dipole moment (μ), quadrupole moment (QT), static dielectric constant (ε0), diffusion constant at infinite dilution (D0,W), viscosity (η), surface potential at the water−air interface in the air to liquid direction (χ, 1 eV = 23.061 kcal/mol), solvent Galvani potential (ØG), and surface tension (γ). Values in parentheses are standard deviations. The standard deviations for χ and ØG amount to 5 mV and 0.1 kcal/mol, respectively. b Dipole moment of liquid water from X-ray diffraction measurements105 (μ = 1.85 D in the gas phase93). cFrom ref 94. dFrom ref 92. eFrom refs 95 and 96. fFrom ref 97. 1041


1042

−81.1 −79.6 −63.8 −58.9 −67.7

2.6 −9.5 −6.6 0.2 −10.0 −6.7 −8.4 −9.1 −8.9 2.8 2.6 −4.8 1.0 0.1 −4.5 −4.1 −4.6 −4.1 2.6 1.1 0.7

−60.3 −81.7 −81.4 −68.3 −58.9 −65.5 2.4 −0.9

−61.1 −80.6 −80.1 −68.8 −59.7 −66.3 2.4 −0.9

c

SPC/E

2.4 −9.7 −6.7 0.0 −10.3 −6.9 −8.6 −9.3 −9.1 2.5 2.4 −4.8 0.6 −0.2 −4.5 −4.2 −5.2 −4.5 2.4 0.9 0.5

c

SPC

−60.0 −82.5 −82.1 −67.9 −58.6 −65.4 2.5 −1.1

2.6 −9.5 −6.5 0.3 −10.1 −6.6 −8.3 −9.0 −8.8 2.9 2.8 −4.7 1.0 0.2 −4.4 −4.3 −4.4 −3.8 2.7 1.1 0.8

c

d

−58.6 −87.8 −88.4 −65.4 −58.6 −66.0 5.0 −3.0

2.6 −11.7 −4.6 0.4 −12.2 −4.5 −12.7 −11.4 −11.7 2.9 2.8 −3.9 0.8 −0.2 −2.9 −2.9 −5.6 −3.8 2.7 1.7 0.7

SPC/Eb

−61.4 −80.7 −80.3 −69.4 −59.9 −66.6 2.6 −1.2

2.6 −9.5 −6.7 0.2 −10.0 −6.8 −8.3 −9.0 −8.8 2.8 2.6 −4.7 1.0 0.1 −4.5 −4.1 −4.6 −4.1 2.6 1.1 0.7

e

OPC3

−60.5 −80.1 −79.7 −68.5 −59.0 −65.8 2.3 −0.4

2.4 −9.7 −6.8 −0.1 −10.1 −6.8 −8.5 −9.1 −9.0 2.4 2.3 −4.7 0.5 −0.3 −4.5 −4.1 −5.2 −4.6 2.4 0.8 0.5

c

TIP3P c

TIP4P

Neutral Side Chains 2.6 2.5 −9.6 −9.9 −6.7 −6.5 0.2 0.1 −10.4 −10.3 −6.8 −6.7 −8.2 −8.6 −9.0 −9.1 −8.8 −9.0 2.7 2.7 2.5 2.6 −5.4 −5.0 0.9 0.8 −0.1 0.1 −4.6 −4.6 −4.4 −4.2 −4.8 −4.6 −4.2 −4.0 2.6 2.6 1.1 1.1 0.6 0.6 Charged Side Chains −59.7 −57.1 −83.0 −85.3 −82.7 −85.0 −67.2 −64.6 −58.5 −56.0 −65.3 −62.7 2.5 4.0 −1.1 −0.5

f

TIP3P-FB

−57.4 −85.7 −85.2 −64.8 −55.9 −62.8 4.1 −0.6

2.6 −9.7 −6.6 0.2 −10.3 −6.7 −8.4 −9.1 −9.0 2.8 2.6 −5.0 1.0 0.1 −4.5 −4.2 −4.4 −4.0 2.6 1.1 0.7

c

TIP4P-Ew

−56.4 −87.7 −87.4 −63.9 −55.0 −61.9 5.5 −0.9

2.5 −10.1 −6.7 0.1 −10.6 −7.0 −8.8 −9.4 −9.2 2.7 2.5 −5.4 0.9 −0.1 −4.6 −4.4 −4.8 −4.2 2.5 1.0 0.5

c

−57.0 −85.1 −88.5 −65.4 −54.5 −61.9 5.5 −0.8

2.6 −10.0 −6.8 0.2 −10.5 −6.9 −8.8 −9.3 −9.2 2.7 2.6 −5.4 0.8 0.0 −4.6 −4.4 −4.9 −4.1 2.5 1.0 0.5

g

−57.9 −84.7 −88.3 −65.6 −55.6 −62.8 5.0 −1.2

2.2 −10.5 −7.1 −0.4 −11.2 −7.6 −9.7 −10.3 −10.2 1.6 1.5 −6.3 −0.3 −1.4 −4.8 −4.9 −6.5 −5.5 1.7 0.9 −0.3

h

TIP4P/2005

−53.6 −80.7 −91.2 −61.6 −53.9 −59.5 6.8 0.9

2.6 −7.3 −5.8 0.1 −11.4 −6.1 −12.5 −5.1 −6.6 2.6 2.6 −3.1 0.5 0.8 −4.8 −3.6 −3.2 −1.5 2.5 2.0 1.3

i

−55.0 −80.3 −91.2 −61.9 −54.9 −60.4 6.4 0.5

2.2 −7.9 −6.3 −0.4 −12.1 −6.8 −13.5 −6.2 −7.6 1.5 1.5 −4.1 −0.7 −0.5 −4.9 −4.1 −4.9 −3.1 1.5 1.4 0.4

j

−65.2 −73.2 −72.8 −75.5 −61.5 −69.4 7.2 −0.2

2.3 −10.2 −8.7 −0.1 −10.5 −8.3 −8.4 −8.8 −8.7 2.4 2.2 −3.8 0.6 0.0 −5.5 −4.9 −5.7 −6.0 2.3 1.0 0.1

c

TIP5P

−65.6 −73.7 −73.1 −75.8 −61.9 −69.5 7.1 −0.6

2.3 −10.3 −8.9 −0.2 −10.7 −8.5 −8.7 −9.0 −9.0 2.2 2.2 −4.0 0.6 −0.2 −5.6 −5.0 −6.1 −6.1 2.2 1.1 −0.1

c

TIP5P-Ew

b

Standard deviations for neutral and charged analogues amount to approximately 0.05 and 0.14 kcal/mol, respectively. RMSE means root-mean-square error, and MSE means mean signed error. Experimental observations for neutral analogues were taken from refs 110 and 111 and those for charged ones from refs 11, 83, and 111 on the basis of ΔGhyd * (H+) = −262.4 kcal/mol. cff99SB-ILDN. d ff15ipq. eff14SB. ffb15. gff99SBw. hff99SBws. iff03w. jff03ws. kA weighted value computed as 0.2·Hid +0.8·Hie.83,112 lGuanidinium (Gud), a partial analogue of Arg+.

a

Arg Asp Glu Gudl His Lys RMSE MSE

1.94 −9.68 −6.7 −1.24 −9.38 −6.47

Ala Asn Asp Cys Gln Glu Hid Hie Hisk Ile Leu Lys Met Phe Ser Thr Trp Tyr Val RMSE MSE

−10.27 2.15 2.28 −4.38 −1.48 −0.76 −5.06 −4.88 −5.88 −6.11 1.99

expt

AA

b

Table 2. Real Hydration Free Energies (in kcal/mol) of Amino Acid Side Chain Analoguesa

Journal of Chemical Information and Modeling Article


Article

Journal of Chemical Information and Modeling water (μ = 2.9 D),105 all of the water models have a smaller μ (by design, since these are empirical models), and no obvious differences are observed among the water models except for TIP4P, which has the smallest value (μ = 2.178 D) (Table 1). The “FB” and “OPC” water models reproduce the experimental dielectric constants (ε0) of water accurately, in perfect agreement with the original reports.28−30 TIP3P, TIP5P, and TIP5P-Ew exhibit larger ε0 than the experiment, whereas the other models produce an underestimation, with TIP4P having the lowest ε0. It should be noted that the ε0 prediction of water models depends highly on the simulation protocol,106 and there is a large diversity in the literature. For the solvent Galvani potential (ØG), which needs to be considered in the calculation of physical (real) solvation free energies (SFEs) for molecular ions, the three- and four-site water models give ØG values ranging from −11.2 to −14.4 kcal/mol/e. This leads to deviations in the calculated ionic SFEs in water from typical MD/PBC simulations by 11−14 kcal/mol (underestimation for monovalent cations and overestimation for monovalent anions). The five-site water models TIP5P and TIP5P-Ew have a small contribution from the surface potential (ØG ≈ −2 kcal/mol/e), which could lead to the conclusion that the water−vacuum interface is not very structured. These models also have a low surface tension, although there is no correlation between the surface potential and the surface tension in the models tested. The predictions of the surface potential (χ) of SPC, SPC/E, and TIP3P in this work are close to the available reports,79,107 while χ for TIP4P (−0.523 V) is less negative than that in the literature (χ ≈ −0.6 V),108,109 leading to a difference of ∼2 kcal/mol/e in ØG. 3.2. Hydration of Amino Acid Side Chains. The hydration free energies (ΔGhyd) of the amino acid side chain analogues were corrected by adding a term based on the solvent’s Galvani potential (qØG in eq 3), and the corrected real HFEs are given in Table 2. The correction amounts to zero for neutral analogues (q = 0). For the neutral amino acid side chain analogues with ff99SB-ILDN, no significant differences are observed among three- and four-site water models, with the root-mean-square error (RMSE) from experimental observations amounting to ∼1 kcal/mol, in line with our previous work for Amber ff99SB, CHARMM 27, GROMOS 54A7/A8, and OPLS-AA/L.83 A similar RMSE (∼1 kcal/mol) for the five-site water models TIP5P and TIP5P-Ew is observed (Table 2); however, these two water models yield more negative ΔGhyd for Thr (T), Trp (W), and Tyr (Y) than the three- and four-site water models by 1−2 kcal/mol, resulting in a good match with the experiments (Table 2 and Figure 3a). The TIP5P and TIP5P-Ew models overestimate the hydration (more negative ΔGhyd) of neutral Asp and Glu residues by ∼2 kcal/mol, whereas the three- and four-site water models reproduce the experimental ΔGhyd values of these two residues accurately. For the three- and four-site water models, in general, ff99SBILDN appears to have a systematic bias in that the hydration of neutral side chains is on average less favorable than experimental values, similar to that with ff99SBw/TIP4P/ 2005, ff14SB/TIP3P, and ff15ipq/SPC/Eb, as indicated by positive mean signed errors (MSEs) of 0.5−0.8 kcal/mol (Tables 2 and S1). These results are in agreement with the work on a partial set of neutral side chain analogues (i.e, most neutral forms of titratable AAs such as Asp, Glu, His, and Lys were not examined) with Amber ff99 and ff03*.17,20 However, TIP5P and TIP5P-Ew show a reasonable representation of amino acid hydration with small MSEs of 0.1 and −0.1 kcal/

Figure 3. Comparison of calculated real hydration free energies using ff99SB-ILDN for (a) neutral and (b) charged amino acid side chain analogues in different water models with experimental observations. Amino acids with large discrepancies in the calculated hydration free energies with varied water models are marked by their one-letter codes (R stands for Gud+ here).

mol, respectively (Table 2). The latest force field, ff15ipq, underestimates ΔGhyd of Asp, Glu, Met, Ser, Thr, and Tyr by about 2 kcal/mol and overestimates those of Asn, Gln, and His by 1−3 kcal/mol, resulting in a large RMSE of 1.7 kcal/mol from experiment. It should be noted that Amber ff99SB-ILDN is identical to ff99SBw for neutral side chain analogues. In order to validate the Amber ff99SBws, HFEs for all of the side chain analogues were calculated with (ff99SBws) and without (ff99SBw) scaling of the protein−water interactions (Table 2). The scaling in ff99SBws helps to strengthen the hydration of the neutral analogues by lowering ΔGhyd for all residues, yielding an overall more favorable hydration with an MSE of −0.3 kcal/mol. Compared with ff99SBw and the experimental observations, ff03w (with parameters for neutral residues identical to those of ff03) displays worse performance with a larger RMSE (∼2 kcal/ mol) and a more positive MSE (1.3 kcal/mol). For instance, the HFEs of Asn, His, Lys, Trp, and Tyr are underestimated by 2−4 kcal/mol in ff03w, as indicated by the neutral analogues (Tables 2 and S2). The scaling of the water interactions in ff03ws produces a more negative ΔGhyd and leads to an overall improvement in RMSE (1.4 kcal/mol) and MSE (0.4 kcal/ mol) relative to ff03w. However, the performance of ff03ws for predicting the hydration of neutral amino acids is still slightly worse than those of ff99SB-ILDN and ff99SBw(s). The intrinsic HFEs (ΔGintr, i.e., without the air−liquid surface contribution in eq 3) of charged amino acids depend critically on the type of water model,83,113 and the difference in ΔGintr among different water models is even larger than 20 kcal/mol. For Glu−, for instance, SPC/E, TIP3P, TIP4P/2005, and TIP5P-Ew with ff99SB-ILDN predict ΔGintr values of −95.0, −92.0, −98.6, and −75.5 kcal/mol, respectively. With the surface contribution, the differences in the real HFEs (ΔGhyd) resulting from the water models are reduced significantly, as indicated by the values for the charged 1043


42.6 33.0 94.7 108.5 108.5

43.5 34.5 96.4 110.8 113.0

0.290 0.290 0.288 0.286 0.288 0.286 0.264 0.264

b 0.302 0.306 0.290 0.290 0.288 0.290 0.262 0.264

c

44.5 36.9 99.3 115.3 118.8

SPC/Eb

43.4 35.5 92.2 104.6 106.8

0.292 0.286 0.288 0.286 0.284 0.284 0.266 0.266

d

OPC3

43.5 34.5 104.1 100.3 104.6

0.290 0.288 0.284 0.282 0.286 0.280 0.268 0.266

b

TIP3P b

TIP4P

RRDF (nm) 0.292 0.294 0.290 0.294 0.288 0.286 0.290 0.288 0.288 0.286 0.286 0.286 0.264 0.264 0.266 0.262 ρ (nm−3) 42.9 41.9 35.5 33.6 97.4 93.4 107.3 109.6 111.9 110.1

e

TIP3P-FB

42.0 35.5 99.8 111.3 110.7

0.292 0.292 0.286 0.288 0.288 0.286 0.264 0.262

b

TIP4P-Ew

41.9 34.1 99.3 113.0 114.8

0.294 0.294 0.292 0.290 0.290 0.290 0.264 0.264

b

41.9 35.5 95.7 100.7 114.7

0.294 0.294 0.292 0.292 0.292 0.288 0.266 0.266

f

TIP4P/2005

40.2 33.6 89.8 96.1 117.0

0.300 0.300 0.296 0.292 0.294 0.290 0.270 0.266

g

57.2 46.0 116.1 88.5 85.9

0.278 0.278 0.276 0.274 0.282 0.278 0.274 0.274

b

TIP5P

58.7 44.9 124.6 89.6 88.7

0.278 0.278 0.276 0.276 0.278 0.278 0.274 0.272

b

TIP5P-Ew

0.285j 0.285k

0.273i

∼0.3h

expt

a Positions of the first peak (RRDF) in the radial distribution functions (RDFs) of water (oxygen atoms) around the ions (nitrogen atoms for cations and oxygen for anions) and water number densities (ρ) at RRDF are given. For pure water, ρ ≈ 33.4 nm−3, corresponding to a density of 1 g/mL. Guandinium (Gud+), imidazolium (Imd+), and ammonium (NH4+) are partial analogues of Arg+, His+, and Lys+ side chains, respectively. bff99SB-ILDN. cff15ipq. dff14SB. efb15. fff99SBws. gff03ws. hEstimate based on neutron diffraction experiments.87 iAverage hydrogen-bond length with water observed in a solid crystal.88 jFrom refs 89 and 90. kFrom X-ray and NMR experiments.91

41.6 33.1 94.4 103.5 102.9

0.290 0.290 0.288 0.290 0.286 0.284 0.266 0.266

Arg+ Gud+ His+ Imd+ Lys+ NH4+ Asp− Glu−

Arg+ His+ Lys+ Asp− Glu−

b

b

AA

0.290 0.286 0.284 0.284 0.284 0.284 0.266 0.264

SPC/E

SPC

Table 3. Characterization of the First Peak of the Ion−Water Radial Distribution Functiona

Journal of Chemical Information and Modeling Article

1044


Article


Compared with the “native” water model of the Amber ff99SB-ILDN force field (TIP3P), TIP4P/2005 slightly overestimates RRDF for water and positively charged amino acid side chains by 0.004−0.008 nm and underestimates RRDF for water and negatively charged amino acid side chains by 0.002−0.004 nm. In contrast, TIP5P and TIP5P-Ew yield large underestimations (which are much further away from the experiments) for positively charged amino acids by 0.006−0.012 nm and overestimations for negatively charged amino acids by ∼0.007 nm (Table 3). RDF plots of Arg+ and Asp−, for instance, are given in Figure 4, pinpointing the discrepancies in

analogues in Table 2. Interestingly, the ΔGhyd values with ff99SB-ILD for the nine examined water models can be putatively classified into three groups on the basis of the number of interacting sites in the water models. The magnitude of ΔGhyd with three-site water models lies between those of four-site water models and five-site water models (Figure 3b). The four-site models overestimate the hydration of negatively charged amino acid side chains but underestimate the hydration of the positively charged side chains; the opposite behavior is observed for the five-site water models. For each group, the difference in ΔGhyd with varied water models for the same amino acid side chain is less than 2 kcal/mol. If the solvent Galvani potential for TIP4P from the literature (ØG ≈ −14 kcal/mol/e)108,109 were to be used instead of the value of −12.1 kcal/mol/e in Table 1, the distinctions between TIP4P and the three-site water models would be not as obvious as in Figure 3b. The debated experimental reference for ΔGhyd of charged residues,83 however, precludes any clear conclusion about the performance of the various water models. The atomic charges of charged amino acids in ff99SB-ILDN, ff99SBw(s), and ff03w(s) are different and thus cannot be compared with each other directly. For the latter two force fields, it should be noted that there is a large difference in the hydration of Asp− and Glu− residues, whereas the difference between the two residues is tiny in the experiments and in the calculations with ff99SB-ILDN and ff15ipq (Table 2). Refitted atomic charges for charged residues in ff99SBw20,36 give accuracy for ΔGhyd comparable to that in ff99SB-ILDN, except for an underestimation of Asp by ∼3 kcal/mol. Scaling of the protein−water interactions in ff99SBws and ff03ws gives a slightly strengthened HFE, and the real HFEs of charged side chains using ff99SBw(s) and ff03w(s) with TIP4P/2005 appear to fall within the group of four-site water models as classified in Figure 3b. The effects of the water model on the hydration of charged amino acids were further investigated in terms of structural ion−water interaction properties. Table 3 lists the locations of the first peaks (RRDF) in the ion−water radial distribution functions (RDFs) and the water number density (ρ) at RRDF, along with, if available, the corresponding experimental estimates. A quick look suggests that the three- and four-site water models reproduce the experimental RRDF values for Gud+ and NH4+ ions well, whereas the five-site water models give an underestimation. The five-site water models perfectly reproduce the experimental RRDF of Imd+, while an overestimation by 0.01 nm is found for the three- and four-site water models. All of the tested water models underestimate RRDF of Asp− by 0.01−0.02 nm. A closer inspection of the relationship between water parameters (Table S3) and water-related properties (Tables 1 and 2) indicates that RRDF has a stronger correlation with the LJ repulsive term (C12, R = 0.6−0.9) of the water model than the dispersive term (C6, R = 0.2−0.6); for the σ−ε form of the LJ potential, RRDF correlates to σ more strongly (R > 0.86) than ε (R < 0.5), as given in Table S4 and Figure S1. RRDF does not show any correlation with the water dipole moment, whereas a strong correlation with the quadrupole moment (QT) and the hydration free energies of charged amino acid side chains is observed. Stronger hydration means a smaller RRDF and a higher ρ at RRDF, as shown in Figure S1 and Table 3. It is found that a reduction in C12, σ, or QT could yield a smaller and larger RRDF for positively and negatively charged side chains, respectively (Figure S1).

Figure 4. Radial distribution functions (RDFs) for (a) nitrogen of Arg+ and (b) oxygen atoms of Asp− around water (HOH) oxygen atoms computed from 10 ns MD simulations with ff99SB-ILDN in TIP3P (solid lines), TIP4P/2005 (dotted), and TIP5P-Ew (dashdotted) water. Unnormalized RDFs are presented here describing water number density (ρ) around charged amino acids; ρ approaches 33.4 nm−3, corresponding to the pure water density of 1 g/mL. Insers are enlarged views of the water number density curves at r = 0.25−0.35 nm.

water−ion binding structures with different water models, which indicates that the five-site models likely produce different hydrated structures from three- and four-site models. With different atomic charges and scaling of the protein− water interactions, ff99SBws and ff03ws with TIP4P/2005 produce results comparable to those with ff99SB-ILDN (Table 3). No significant differences for the Amber FFs ff14SB, ff15ipq, and fb15 were observed, except that ff15ipq overestimates the RRDF of Gud+ by 0.006 nm. 3.3. Diffusion of Amino Acids. Similar to water diffusion constant calculations, translational diffusion coefficients of zwitterionic amino acids at infinite dilution were determined by extrapolation according to eq 2 (see Table 4). With ff99SBILDN, TIP3P yields the largest RMSE from experimental observations (1.44 × 10−5 cm2/s), while SPC/Eb gives the smallest RMSE (0.15 × 10−5 cm2/s), despite the fact that SPC/ Eb underestimates the diffusion constant by ∼9% (Table 1). The RMSEs with various water models are found to be correlated strongly (R = 0.99) with the water self-diffusion constants (D0), as shown in Figure 2; the faster the water 1045


Article


Table 4. Comparison of Calculated Translational Diffusion Constants (in 10−5 cm2/s) of Amino Acids with Experimental Observations at Infinite Dilutiona AA

exptb

Ala 0.74 Arg 0.55 Asn 0.66 Asp 0.65 Cys 0.69 Gln Glu 0.62 Gly His 0.56 Ile 0.62 Leu Lys 0.53 Met 0.64 Phe 0.60 Ser 0.71 Thr Trp Tyr Val RMSEcalci RMSEscaledi R

SPC

SPC/E

OPC3

TIP3P

TIP3P-FB

TIP4P

TIP4P-Ew

c

c

c

SPC/Eb d

e

c

f

c

c

c

TIP4P/2005 g

h

TIP5P c

TIP5P-Ew c

1.87 1.63 1.56 1.52 1.71 1.52 1.46 2.17 1.38 1.52 1.65 1.56 1.54 1.40 1.68 1.66 1.59 1.43 1.73 0.94 0.18 0.66

1.09 0.99 0.87 1.16 1.13 1.05 0.91 1.44 1.13 0.97 1.04 0.95 1.16 0.94 1.04 1.32 0.84 0.99 0.99 0.41 0.24 0.29

0.79 0.74 0.82 0.63 0.92 0.69 0.69 0.92 0.74 0.78 0.78 0.71 0.75 0.80 0.78 0.79 0.74 0.78 0.83 0.15 0.21 0.41

0.75 0.66 0.77 0.69 0.85 0.83 0.70 0.93 0.75 0.84 0.72 0.75 0.78 0.83 0.92 0.76 0.70 0.64 0.81 0.16 0.23 0.44

0.99 0.72 0.86 0.78 0.92 0.84 0.81 1.05 0.87 0.84 0.80 0.84 0.90 0.92 0.94 0.89 0.83 0.88 0.82 0.24 0.20 0.65

2.40 1.96 2.16 1.89 2.33 1.98 1.83 2.76 2.00 2.07 2.17 1.75 2.22 1.89 2.20 2.30 1.83 1.82 1.87 1.44 0.14 0.80

0.87 0.73 0.87 0.70 0.76 0.78 0.71 1.13 0.87 0.71 0.84 0.74 0.78 0.83 0.95 0.75 0.67 0.74 0.72 0.18 0.22 0.43

1.70 1.25 1.32 1.43 1.51 1.39 1.39 1.82 1.48 1.34 1.47 1.48 1.53 1.30 1.56 1.41 1.31 1.32 1.35 0.82 0.24 0.59

1.07 0.97 0.95 0.93 1.00 0.92 1.10 1.36 1.07 0.95 0.99 0.91 1.01 0.88 1.07 1.07 0.83 0.84 0.93 0.37 0.25 0.39

0.96 0.82 0.88 0.84 0.93 0.89 0.72 1.01 0.86 0.92 0.78 0.88 0.87 0.80 1.01 0.95 0.79 0.86 0.85 0.25 0.25 0.57

0.93 0.83 0.81 0.74 1.04 0.91 0.75 1.16 0.84 0.78 0.78 0.91 0.96 0.79 1.02 0.84 0.78 0.81 0.81 0.25 0.25 0.43

0.98 0.80 0.81 0.80 0.99 0.68 0.67 1.10 0.86 0.78 0.87 0.76 0.75 0.88 0.80 0.85 0.81 0.82 0.90 0.21 0.21 0.45

0.97 0.87 0.93 0.95 0.98 0.78 0.87 1.16 0.96 1.08 0.75 0.87 0.93 1.06 0.99 0.64 0.96 0.85 1.01 0.33 0.13 0.33

1.24 0.96 0.94 0.94 1.00 1.10 0.92 1.23 1.05 1.01 1.00 0.88 1.03 0.88 0.95 0.84 0.85 0.94 1.01 0.36 0.12 0.52

Amino acids are in zwitterionic forms for consistency with experimental determinations at pH ∼3.5; relative errors for the calculations are ∼7%. Taken from ref 52. Standard deviations for the experiments are smaller than 0.01 × 10−5 cm2/s and are not given here. cff99SB-ILDN. dff15ipq. e scaled ff14SB. ffb15. gff99SBws. hff03ws. iRoot-mean-square errors (RMSEs) with respect to experiment of the calculated (Dcalc 0,AA) and scaled (D0,AA ; eq 4) diffusion constants. a b

(Figure 5). Because of the similar results for water diffusivity, ff15ipq/SPC/Eb performs similarly to ff99SB-ILDN/SPC/Eb and fb15/TIP3P-FB, and ff14SB/OPC3 performs similarly to ff99SB-ILDN/TIP4P/2005 (Table 4). It should be noted that ff99SB-ILDN and ff99SBw have identical force field parameters for zwitterionic amino acids, as stated in the model construction in section 2.1. Unexpectedly, scaling of the water interactions in ff99SBws has little influence on the diffusion properties of amino acids (Table 4), although such scaling leads to more favorable hydration of amino acid side chains (Table 2). Also, Amber ff03ws yields diffusion results comparable to those with ff99SB-ILDN and ff99SBws calc exp (Table 4). With a scaling factor of Dexp 0,W/D0,W, where D0,W and calc D0,W are the experimental and calculated water diffusion coefficients (eq 4), comparable results for all of the water models tested can be obtained.

model diffuses, the larger the RMSE is. The relationship of the calculated and experimental diffusion constants of the amino acids with their molar masses is given in Figure 5. All of the

■

Figure 5. Calculated and experimental diffusion constants of amino acids at infinite dilution as functions of amino acid molar mass with ff99SB-ILDN. Solid lines are fits to a power function indicative of a relationship with the molar mass. Data points for TIP4P-Ew, TIP5P, and TIP5P-Ew are close to those of SPC/E and TIP4P/2005 and for clarity are not given here.

DISCUSSION A systematic benchmark of protein force fields (FFs) for reproducing hydration properties and diffusion of amino acids is presented to assess the performance of FF/water model combinations. Eight Amber force fields, namely, ff99SBILDN,10 ff99SBw,20 ff99SBws,20 ff03w,33 ff03ws,20 ff14SB,50 ff15ipq,18 and fb15,41 were examined with 11 water models belonging to three categories, namely, three-site (OPC3,30 SPC,14 SPC/E,15 SPC/Eb,19 TIP3P,13 and TIP3P-FB28), foursite (TIP4P,13 TIP4P-Ew,25 and TIP4P/200527), and five-site (TIP5P16 and TIP5P-Ew26) models. The widely used variants of Amber-type fixed-charge force fields, as examined in this work, share the potential form and most parameters from the original ff94 model by Cornell et al.,114 where atomic charges

water models examined overestimate the diffusion constants of zwitterionic amino acids with ff99SB-ILDN. The fast diffusion of TIP3P greatly accelerates the diffusion of amino acids and produces the overall largest diffusion constants, which are 2 times larger than the experimental observations. TIP4P/2005, which performs best for the water diffusivity, still gives an RMSE of 0.25 × 10 −5 cm 2 /s, an overestimation by approximately 40%. As expected, diffusion of amino acids decreases, probably in a power way, as the molar mass increases 1046


Article


Figure 6. Comparison of (a, b) dipole moments, (c, d) atomic charges, and (e, f) hydration free energies in TIP3P for the constructed neutral amino acid side chain analogues by the method of Cβ modification and QM full ESP calculations for (top) Amber ff99SB-ILDN and (bottom) ff03.

were obtained by fitting the calculated electrostatic potential (ESP) at the HF/6-31G* level of theory in the gas phase with the restrained electrostatic potential (RESP) method.115 Subsequent adaptations were made, focusing primarily on the improved representation of protein secondary structures by successive tuning of the ff94 dihedral torsions (e.g., ff99116 and ff99SB35) and side chain torsions (ff99SB-ILDN10). In the latest Amber FF, ff14SB, for instance, all of the amino acid side chain dihedral parameters that were carried over from ff94 were completely refitted,50 while Amber fb15 started from ff99SB, modified the bonded parameters by fitting to high-quality quantum-mechanical (QM) calculations with the force balance (FB) method,28 and left the other parameters unmodified.41 Inheriting from the same parent (ff94), ff99SB-ILDN, ff14SB, and fb15 therefore share identical atomic charges and nonbonded parameters for all amino acids, of which charged residues have the same backbone charges and the neutral ones share a different set of backbone charges. In order to eliminate the difference in backbone charges for a faithful reproduction of helix propensities, Best and Hummer53 refitted the atomic charges of charged residues in ff99SB-ILDN and side chain torsions for Asp in conjunction with the new partial charges, and together with a fine-tuning of backbone dihedral potentials, the resulting force field is referred to as ff99SB*-ILDN-Q.36 On this basis, ff99SBws was developed in order to get strengthened protein−water interactions.20 HF/6-31G* has been reported to overestimate the gas-phase dipole moment and to produce slightly larger charges, leading to an overestimated gas-phase dipole moment, which is somewhat condensed-phase-like.54 For a more accurate representation of polarization effects, the Amber ff03 lineage applied the ESP for RESP fitting with an increased QM level of theory (B3LYP/cc-pVTZ) using the IEFPCM continuum solvent model (with an effective dielectric constant of 4, mimicking organic media or protein interiors).54 Unlike ff94, each amino acid in ff03, on which ff03ws was based, has its own backbone charges, although bonded (except main-chain torsion) and van der Waals parameters were still taken from ff94 directly. The third charge model in the evolution of Amber protein models was the implicitly polarized charge (IPolQ) scheme,117 which assigns polarized atomic charges implicitly in

the presence of explicit solvent, as was done in the latest Amber FF, ff15ipq, with the SPC/Eb water model.18 ff15ipq is an original protein FF with its own charge model and new angle and torsion terms as well as new nonbonded parameters for the polar hydrogens designed to improve amino acid salt bridge associations.18 The consistency used when treating backbone charges on the peptide group automatically compartmentalizes the charge representations and thereby lends credence to a manual charge assignment. The procedure introduced by Shirts et al.48 for taping together representations of amino acid side chain analogues was therefore used in this work as well and is called the method of Cβ modification. However, ff03 does not have any true consistency in backbone charges, which makes manual charge assignment somewhat arbitrary. Following the original papers,54,114 a full ESP charge reassignment by QM calculations (called the method of QM full ESP) was attempted in order to construct the side chain analogues and zwitterionic amino acids with the Gaussian 09 software.118 Compared with QM full ESP, Cβ modification slightly overestimates the dipole moment of neutral side chain analogues (Figure 6a,b) and underestimates the corresponding atomic charges (in absolute value) by ∼20% (Figure 6c,d). Interestingly, the QM full ESP method does not seem to influence the hydration free energies (ΔGhyd) of the analogues significantly for ff99SB-ILDN (Figure 6e and Table S5), but it shows a better prediction for ff03 than the manual assignment by Cβ modification, as indicated by a reduced MSE from 1.3 to 0.8 kcal/mol and a reduced RMSE from 1.8 to 1.0 kcal/mol (Figure 6f and Table S2), a similar performance to that with f99SB-ILDN. That is, the ΔGhyd performance of ff03 is better than that observed in Table 2. Discrepancies between the manual assignment and QM full ESP calculations were observed as well for the dipole moments and atomic charges of zwitterionic amino acids (Figure S2). These findings indicate that the charge distribution could be a source of error for such calculations. Our benchmark results for neutral amino acid side chains and neutral zwitterionic amino acids are summarized in Figure 7. The “native” model TIP3P in combination with ff99SB-ILDN shows the smallest RMSE of 0.8 kcal/mol from experiment for the hydration free energies of the neutral analogues. The RMSE 1047


Article


molecules should be developed and optimized for use with the chosen water model.119−121 Despite the good performance on ΔGhyd, TIP3P greatly overestimates amino acid diffusion constants, which is likely due to the inherent deficiency of the TIP3P water model, in particular, the significantly overestimated self-diffusivity (D0 = (6.14 ± 0.05) × 10−5 cm2/s) and the corresponding low viscosity (η = 0.32 ± 0.04 mPa·s) compared with the experimental observations (2.299 × 10−5 cm2/s and 0.89 mPa·s, respectively). With the optimized TIP4P/2005 model, which accurately captures many properties of liquid water, including diffusion constant (D0 = 2.31 ± 0.02 × 10−5 cm2/s) and viscosity (η = 0.84 ± 0.08 mPa·s), better performance with ff99SB-ILDN (i.e., much closer to the experiments) for absolute diffusion constants of amino acids is obtained, similar to the predictions by the latest OPC3 water model with Amber ff14SB (Table 4). Using a water model with a lower D0 than experiment, like SPC/Eb, reduces the diffusivity of amino acids further, but a tendency of overestimation still exists for TIP4P/ 2005 and SPC/Eb as well as the other water models tested, even with the most recent FFs, ff14SB, ff15ipq, and fb15 (Table 4 and Figure 5). With an ad hoc correction for the diffusion properties of the water model (eq 4), the calculated results for different water models are comparable (Table 4). These findings indicate an improvement in the simulated diffusivity of amino acids with an accurate model for liquid water and the necessity of improving the accuracy of the force fields themselves. In view of the fact that despite the very accurate diffusion properties of TIP4P/2005, the diffusion constants of amino acids are still overestimated by about 40%, even when using the scaled water interactions in ff99SBws and ff03ws,20 something still seems to be missing from the force fields. A modification of charges in solvated side chains (large absolute values of the charges) would probably reduce the diffusion rates, but most likely at the cost of overestimated HFEs. It has been suggested, on the basis of extensive simulations of molecular liquids, that the dispersion interactions in biomolecular force fields are overestimated104 because the force fields are optimized without long-range dispersion corrections. Although a reduction of the dispersion constant (C6 in the Lennard-Jones potential) could have an effect on protein properties, it is unlikely to affect the diffusion of isolated amino acids. It has been suggested that the balance between protein− water interactions and protein−protein interactions for current biomolecular force fields is not correct for the modeling of disordered proteins and protein aggregates.20,42−46 The slightly too weak protein−water interactions are confirmed here by positive MSEs for the hydration of neutral side chain analogues with ff99SB-ILDN, ff99w, ff03w, ff14SB, ff15ipq, and fb15 (Table 2 and Figure 7). Tuning the solute−water van der Waals interactions appears to be a partial solution to this issue,20,46,49 as was done in the ff99SBws and ff03ws force fields by Best et al.20 For these two force fields, scaling of the protein−water interactions was proposed to modulate (strengthen, more precisely) the protein−water interactions and generate, likely without highly residue-specific effects, more negative hydration free energies of amino acids, as confirmed by reduced MSEs from 0.5 to −0.3 kcal/mol for ΔGhyd of neutral side chains with ff99SBws and from 1.3 to 0.4 kcal/mol with ff03ws (Table 2 and Figure 7). More importantly, such scaling shows improved properties of simulated disordered proteins and nonspecific protein association.20 The strengthened interaction between

Figure 7. Overall performance of the examined force field/water combinations as indicated by (a) root-mean-square error (RMSE) and (b) mean signed error (MSE) for hydration free energies (ΔGhyd) of neutral amino acid side chains as well as by the MSE for diffusion constants (D0,AA) of neutral zwitterionic amino acids at infinite dilution compared with experiment. Because all of the water models overestimate amino acid diffusion constants, the RMSEs for D0,AA are almost identical to the corresponding MSEs.

amounts to about 1 kcal/mol for all the FFs that were inherited from ff94 and shared an identical charge model, whereas ff03w and ff15ipq yield large RMSEs of 1.4 and 1.7 kcal/mol, respectively (Figure 7a). It should be noted that the QM full ESP charge reassignment leads to a reduced RMSE for ff03w, as mentioned above in Figure 6f. ΔGhyd of the neutral analogues for ff99SB-ILDN paired with three- and four-site models is almost independent of the water model used, with MSEs of 0.5−0.8 kcal/mol, which are less favorable compared with experiment. However, TIP5P and TIP5P-Ew tend to give a more negative ΔGhyd with ff99SB-ILDN for certain amino acid residues, which is better for some but worse for others compared with the other water models and experiments, resulting in overall small MSEs of 0.1 and −0.1 kcal/mol, respectively (Table 2 and Figures 3a and 7b). Interestingly, ΔGhyd for the charged analogues with ff99SBILDN is water-model-dependent and can be classified putatively into three groups, depending on the number of interacting sites in the water model (Figure 3b), of which the three-site models perform best with an RMSE of 2 kcal/mol. It should be noted that such a classification is based on the comparison with experiments with a reference ΔG*hyd(H+) value of −262.4 kcal/mol.85 For each group, the difference in ΔGhyd with different water models for the same amino acid side chain is less than 2 kcal/mol and appears to be negligible, and no obvious differences in the locations of first hydration peak (RRDF) in the radial distribution functions between charged amino acids and water molecules are observed (Table 3). Even so, however, the choice of water model was reported to show a strong influence on the accuracy of other properties of amino acid side chains such as hydration entropies, enthalpies, and heat capacities.17 If other molecules such as metal ions are involved in the simulations, force field parameters for these 1048



■

solute and water molecules is expected to slow down the diffusion of the solute; however, this scaling does not influence the amino acid diffusion significantly (Table 4; see the above discussion). The location of the first hydration peak (RRDF) also needs to be considered for further improvements of the protein force fields. The three- and four-site water models with Amber force fields largely overestimate RRDF of Imd+ by 0.01 nm and underestimate RRDF of Asp− by 0.01−0.02 nm. Except for Imd+, whose RRDF is reproduced accurately, the five-site water model appears to predict a systematic underestimation of RRDF for charged amino acids (Table 3). Such deficiencies likely lead to a deviation in the simulated biomolecular structures from the experiments. In practice, the water−ion distance could probably be fixed by modifying the protein force field, such as the efforts made in the GROMOS 54A8 parameter set,11 or by refinement of the water model. A reduction in the repulsion term (C12), σ in the σ−ε form of the Lennard-Jones potential, or the quadrupole moment (QT) of the water model, for instance, could help somewhat to obtain better hydrated structures (Figure S1), but with necessary adjustments for reproducing key properties of water molecules. The benchmark results highlight the complexity of modeling amino acid hydration with different water models, in particular for charged residues. Atomic charges in the Amber ff94 and ff03 lineages appear to perform better for hydration free energies of side chains than those in the most recent FF, ff15ipq, although the implicitly polarized charge (IPolQ) method in ff15ipq made an important leap forward in terms of solvent-induced polarization for condensed-phase simulations. The exaggerated self-diffusivity of the water models is shown to be largely responsible for the observed discrepancy in the predicted diffusion constants of amino acids by various water models. Because of the accurate diffusion properties, TIP4P/2005 and the “FB” and “OPC” families of the most recent generation are recommended for use with applications where the absolute diffusivity or kinetic properties are of interest. The comparison among different water models not only indicates that the use of accurate water models could result in, to some extent, an improvement of the protein force fields but also provides a source of discrepancies/errors in the simulation studies. The limited success (Tables 1 and 4) of more complex four- and five-site water models (although it should be noted that the latest OPC and TIP4P-FB were not tested here) underscores the difficulty of building an accurate water model that simultaneously is a good solvent model for proteins. With simple atom-based models, parameter development for water and proteins will probably lead to compromises in accuracy. Improvement of accuracy will be limited until fast codes and parameter development pipelines can support more advanced representations. Some promising work on polarizable models is underway.122−124 This work establishes a fundamental assessment of the hydration and diffusion of amino acids and provides implications for further improvements and applications of protein force fields. Reasonable modeling of the surrounding environments (other than water) around biomolecules is of critical importance and needs to be checked carefully before adaptation and reoptimization of the biomolecular force fields. As a follow-up, properties of amino acids in organic solvents other than water will be presented in a subsequent article.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.8b00026. Supplementary tables and figures (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Haiyang Zhang: 0000-0002-2410-7078 Chunhua Yin: 0000-0002-5067-2951 Yang Jiang: 0000-0003-1100-9177 David van der Spoel: 0000-0002-7659-8526 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. Tianwei Tan for a grant of computer time through ChemCloudComputing of Beijing University of Chemical Technology. This work was supported by the National Natural Science Foundation of China (21606016), the Beijing Natural Science Foundation (5174036), and the Fundamental Research Funds for the Central Universities (FRF-TP-17-009A2) as well as by the Swedish Research Council (Grant 2013-5947).

■

REFERENCES

(1) Bellissent-Funel, M.-C.; Hassanali, A.; Havenith, M.; Henchman, R.; Pohl, P.; Sterpone, F.; van der Spoel, D.; Xu, Y.; Garcia, A. E. Water Determines the Structure and Dynamics of Proteins. Chem. Rev. 2016, 116, 7673−7697. (2) Jungwirth, P. Biological Water or Rather Water in Biology? J. Phys. Chem. Lett. 2015, 6, 2449−2451. (3) Finney, J. L. The Water Molecule and Its Interactions: The Interaction between Theory, Modelling, and Experiment. J. Mol. Liq. 2001, 90, 303−312. (4) Guillot, B. A Reappraisal of What We Have Learnt During Three Decades of Computer Simulations on Water. J. Mol. Liq. 2002, 101, 219−260. (5) Jorgensen, W. L.; Tirado-Rives, J. Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6665−6670. (6) Wallqvist, A.; Mountain, R. D. Molecular Models of Water: Derivation and Description. Rev. Comput. Chem. 2007, 13, 183−247. (7) Ouyang, J. F.; Bettens, R. Modelling Water: A Lifetime Enigma. Chimia 2015, 69, 104−111. (8) Cisneros, G. A.; Wikfeldt, K. T.; Ojamäe, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bartok, A. P.; Csanyi, G.; Molinero, V.; Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to ManyBody Potential Energy Functions. Chem. Rev. 2016, 116, 7501−7528. (9) Chaplin, M. www.lsbu.ac.uk/water/water_models.html (accessed April 12, 2018). (10) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (11) Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. New Interaction Parameters for Charged Amino Acid Side Chains in the Gromos Force Field. J. Chem. Theory Comput. 2012, 8, 3705−3723. (12) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins Via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474−6487.

1049


Article


Induced Unfolded State Collapse. J. Phys. Chem. B 2010, 114, 14916−14923. (34) Nerenberg, P. S.; Head-Gordon, T. Optimizing Protein-Solvent Force Fields to Reproduce Intrinsic Conformational Preferences of Model Peptides. J. Chem. Theory Comput. 2011, 7, 1220−1230. (35) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (36) Best, R. B.; de Sancho, D.; Mittal, J. Residue-Specific α-Helix Propensities from Molecular Simulation. Biophys. J. 2012, 102, 1462− 1467. (37) Knott, M.; Best, R. B. A Preformed Binding Interface in the Unbound Ensemble of an Intrinsically Disordered Protein: Evidence from Molecular Simulations. PLoS Comput. Biol. 2012, 8, e1002605. (38) Zerze, G. H.; Mittal, J.; Best, R. B. Diffusive Dynamics of Contact Formation in Disordered Polypeptides. Phys. Rev. Lett. 2016, 116, 068102. (39) Miller, C.; Zerze, G. l. H.; Mittal, J. Molecular Simulations Indicate Marked Differences in the Structure of Amylin Mutants, Correlated with Known Aggregation Propensity. J. Phys. Chem. B 2013, 117, 16066−16075. (40) Mittal, J.; Yoo, T. H.; Georgiou, G.; Truskett, T. M. Structural Ensemble of an Intrinsically Disordered Polypeptide. J. Phys. Chem. B 2013, 117, 118−124. (41) Wang, L.-P.; McKiernan, K. A.; Gomes, J.; Beauchamp, K. A.; Head-Gordon, T.; Rice, J. E.; Swope, W. C.; Martínez, T. J.; Pande, V. S. Building a More Predictive Protein Force Field: A Systematic and Reproducible Route to AMBER-FB15. J. Phys. Chem. B 2017, 121, 4023−4039. (42) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States. J. Phys. Chem. B 2015, 119, 5113−5123. (43) Rauscher, S.; Gapsys, V.; Gajda, M. J.; Zweckstetter, M.; de Groot, B. L.; Grubmüller, H. Structural Ensembles of Intrinsically Disordered Proteins Depend Strongly on Force Field: A Comparison to Experiment. J. Chem. Theory Comput. 2015, 11, 5513−5524. (44) Abriata, L. A.; Dal Peraro, M. Assessing the Potential of Atomistic Molecular Dynamics Simulations to Probe Reversible Protein−Protein Recognition and Binding. Sci. Rep. 2015, 5, 10549. (45) Petrov, D.; Zagrovic, B. Are Current Atomistic Force Fields Accurate Enough to Study Proteins in Crowded Environments? PLoS Comput. Biol. 2014, 10, e1003638. (46) Nawrocki, G.; Wang, P.-h.; Yu, I.; Sugita, Y.; Feig, M. SlowDown in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation. J. Phys. Chem. B 2017, 121, 11072− 11084. (47) Shirts, M. R.; Pande, V. S. Solvation Free Energies of Amino Acid Side Chain Analogs for Common Molecular Mechanics Water Models. J. Chem. Phys. 2005, 122, 134508. (48) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely Precise Free Energy Calculations of Amino Acid Side Chain Analogs: Comparison of Common Molecular Mechanics Force Fields for Proteins. J. Chem. Phys. 2003, 119, 5740−5761. (49) Nerenberg, P. S.; Jo, B.; So, C.; Tripathy, A.; Head-Gordon, T. Optimizing Solute−Water Van Der Waals Interactions to Reproduce Solvation Free Energies. J. Phys. Chem. B 2012, 116, 4524−4534. (50) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (51) Zhang, H.; Lv, Y.; Tan, T.; van der Spoel, D. Atomistic Simulation of Protein Encapsulation in Metal−Organic Frameworks. J. Phys. Chem. B 2016, 120, 477−484. (52) Germann, M. W.; Turner, T.; Allison, S. A. Translational Diffusion Constants of the Amino Acids: Measurement by Nmr and Their Use in Modeling the Transport of Peptides. J. Phys. Chem. A 2007, 111, 1452−1455.

(13) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (14) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; Springer: Dordrecht, The Netherlands, 1981; pp 331−342. (15) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (16) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910− 8922. (17) Hess, B.; van der Vegt, N. F. Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic Comparison of Biomolecular Force Fields and Water Models. J. Phys. Chem. B 2006, 110, 17616−17626. (18) Debiec, K. T.; Cerutti, D. S.; Baker, L. R.; Gronenborn, A. M.; Case, D. A.; Chong, L. T. Further Along the Road Less Traveled: Amber ff15ipq, an Original Protein Force Field Built on a SelfConsistent Physical Model. J. Chem. Theory Comput. 2016, 12, 3926− 3947. (19) Takemura, K.; Kitao, A. Water Model Tuning for Improved Reproduction of Rotational Diffusion and NMR Spectral Density. J. Phys. Chem. B 2012, 116, 6279−6287. (20) Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein−Water Interactions Improve Properties of Disordered Proteins and NonSpecific Protein Association. J. Chem. Theory Comput. 2014, 10, 5113− 5124. (21) Yeh, I.-C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108, 15873−15879. (22) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. A Systematic Study of Water Models for Molecular Simulation: Derivation of Water Models Optimized for Use with a Reaction Field. J. Chem. Phys. 1998, 108, 10220−10230. (23) Panteva, M. T.; Giambaşu, G. M.; York, D. M. Comparison of Structural, Thermodynamic, Kinetic and Mass Transport Properties of Mg2+ Ion Models Commonly Used in Biomolecular Simulations. J. Comput. Chem. 2015, 36, 970−982. (24) Rhee, Y. M.; Sorin, E. J.; Jayachandran, G.; Lindahl, E.; Pande, V. S. Simulations of the Role of Water in the Protein-Folding Mechanism. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 6456−6461. (25) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665−9678. (26) Rick, S. W. A Reoptimization of the Five-Site Water Potential (TIP5P) for Use with Ewald Sums. J. Chem. Phys. 2004, 120, 6085− 6093. (27) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (28) Wang, L.-P.; Martinez, T. J.; Pande, V. S. Building Force Fields: An Automatic, Systematic, and Reproducible Approach. J. Phys. Chem. Lett. 2014, 5, 1885−1891. (29) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5, 3863− 3871. (30) Izadi, S.; Onufriev, A. V. Accuracy Limit of Rigid 3-Point Water Models. J. Chem. Phys. 2016, 145, 074501. (31) Beauchamp, K. A.; Lin, Y.-S.; Das, R.; Pande, V. S. Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse Nmr Measurements. J. Chem. Theory Comput. 2012, 8, 1409−1414. (32) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647. (33) Best, R. B.; Mittal, J. Protein Simulations with an Optimized Water Model: Cooperative Helix Formation and Temperature1050


Article

Journal of Chemical Information and Modeling (53) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix−Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113, 9004−9015. (54) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J; Wang, J.; Kollman, P. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (55) Case, D. A.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Luo, R.; Madej, B.; Mermelstein, D.; Merz, K. M.; Monard, G.; Nguyen, H.; Nguyen, H. T.; Omelyan, A. O.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; Botello-Smith, W. M.; Swails, J.; Walker, R. C.; Wolf, R. M.; Wu, X.; Xiao, L.; Kollman, P. AmberTools 16; University of California: San Francisco, 2016. (56) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (57) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (58) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (59) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. (60) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (61) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (62) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (63) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (64) Nose, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (65) Miyamoto, S.; Kollman, P. A. SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (66) Spångberg, D.; Hermansson, K. Effective Three-Body Potentials for Li+(aq) and Mg2+(aq). J. Chem. Phys. 2003, 119, 7263−7281. (67) Tazi, S.; Boţan, A.; Salanne, M.; Marry, V.; Turq, P.; Rotenberg, B. Diffusion Coefficient and Shear Viscosity of Rigid Water Models. J. Phys.: Condens. Matter 2012, 24, 284117. (68) Dünweg, B.; Kremer, K. Molecular Dynamics Simulation of a Polymer Chain in Solution. J. Chem. Phys. 1993, 99, 6983−6997. (69) Hess, B. Determining the Shear Viscosity of Model Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116, 209−217. (70) Einstein, A. On the Motion of Small Particles Suspended in Liquids at Rest Required by the Molecular-Kinetic Theory of Heat. Ann. Phys. 1905, 322, 549−560. (71) R Core Team. R: A Language and Environment for Statistical Computing, version 3.2.3; R Foundation for Statistical Computing: Vienna, Austria, 2015. (72) Neumann, M.; Steinhauser, O.; Pawley, G. S. Consistent Calculation of the Static and Frequency-Dependent Dielectric Constant in Computer Simulations. Mol. Phys. 1984, 52, 97−113. (73) Neumann, M. Dipole Moment Fluctuation Formulas in Computer Simulations of Polar Systems. Mol. Phys. 1983, 50, 841− 858.

(74) Abascal, J. L. F.; Vega, C. The Water Forcefield: Importance of Dipolar and Quadrupolar Interactions. J. Phys. Chem. C 2007, 111, 15811−15822. (75) van der Spoel, D.; Wensink, E. J.; Hoffmann, A. C. Lifting a Wet Glass from a Table: A Microscopic Picture. Langmuir 2006, 22, 5666− 5672. (76) Van Gunsteren, W.; Berendsen, H. A Leap-Frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1, 173−185. (77) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. Lincs: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (78) Zhang, H.; Jiang, Y.; Yan, H.; Cui, Z.; Yin, C. Comparative Assessment of Computational Methods for Free Energy Calculations of Ionic Hydration. J. Chem. Inf. Model. 2017, 57, 2763−2775. (79) Lamoureux, G.; Roux, B. Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J. Phys. Chem. B 2006, 110, 3308−3322. (80) Misin, M.; Fedorov, M. V.; Palmer, D. S. Hydration Free Energies of Molecular Ions from Theory and Simulation. J. Phys. Chem. B 2016, 120, 975−983. (81) Simonson, T.; Hummer, G.; Roux, B. Equivalence of M- and PSummation in Calculations of Ionic Solvation Free Energies. J. Phys. Chem. A 2017, 121, 1525−1530. (82) Hünenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; Royal Society of Chemistry: Cambridge, U.K., 2011. (83) Zhang, H.; Jiang, Y.; Yan, H.; Yin, C.; Tan, T.; van der Spoel, D. Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton. J. Phys. Chem. Lett. 2017, 8, 2705−2712. (84) Li, P.; Song, L. F.; Merz, K. M., Jr Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645−1657. (85) Randles, J. The Real Hydration Energies of Ions. Trans. Faraday Soc. 1956, 52, 1573−1581. (86) Zhang, H.; Ge, C.; van der Spoel, D.; Feng, W.; Tan, T. Insight into the Structural Deformations of Beta-Cyclodextrin Caused by Alcohol Cosolvents and Guest Molecules. J. Phys. Chem. B 2012, 116, 3880−3889. (87) Mason, P.; Neilson, G.; Dempsey, C.; Barnes, A.; Cruickshank, J. The Hydration Structure of Guanidinium and Thiocyanate Ions: Implications for Protein Stability in Aqueous Solution. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 4557−4561. (88) Steiner, T. The Hydrogen Bond in the Solid State. Angew. Chem., Int. Ed. 2002, 41, 48−76. (89) Jensen, K. P.; Jorgensen, W. L. Halide, Ammonium, and Alkali Metal Ion Parameters for Modeling Aqueous Solutions. J. Chem. Theory Comput. 2006, 2, 1499−1509. (90) Marcus, Y. Ionic Radii in Aqueous Solutions. Chem. Rev. 1988, 88, 1475−1498. (91) Caminiti, R.; Cucca, P.; Monduzzi, M.; Saba, G.; Crisponi, G. Divalent Metal−Acetate Complexes in Concentrated Aqueous Solutions. An X-Ray Diffraction and NMR Spectroscopy Study. J. Chem. Phys. 1984, 81, 543−551. (92) Mills, R. Self-Diffusion in Normal and Heavy Water in the Range 1−45. Deg. J. Phys. Chem. 1973, 77, 685−688. (93) Badyal, Y.; Saboungi, M.-L.; Price, D.; Shastri, S.; Haeffner, D.; Soper, A. Electron Distribution in Water. J. Chem. Phys. 2000, 112, 9206−9208. (94) Pátek, J.; Hrubý, J.; Klomfar, J.; Součková, M.; Harvey, A. H. Reference Correlations for Thermophysical Properties of Liquid Water at 0.1 MPa. J. Phys. Chem. Ref. Data 2009, 38, 21−29. (95) Harris, K. R.; Woolf, L. A. Temperature and Volume Dependence of the Viscosity of Water and Heavy Water at Low Temperatures. J. Chem. Eng. Data 2004, 49, 1064−1069. (96) Harris, K. R.; Woolf, L. A. Temperature and Volume Dependence of the Viscosity of Water and Heavy Water at Low Temperatures. J. Chem. Eng. Data 2004, 49, 1851−1851. 1051


Article

Journal of Chemical Information and Modeling (97) Vargaftik, N.; Volkov, B.; Voljak, L. International Tables of the Surface Tension of Water. J. Phys. Chem. Ref. Data 1983, 12, 817−820. (98) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. Prediction of SelfDiffusion Coefficient and Shear Viscosity of Water and Its Binary Mixtures with Methanol and Ethanol by Molecular Simulation. J. Chem. Phys. 2011, 134, 074508. (99) Mao, Y.; Zhang, Y. Thermal Conductivity, Shear Viscosity and Specific Heat of Rigid Water Models. Chem. Phys. Lett. 2012, 542, 37− 41. (100) González, M. A.; Abascal, J. L. The Shear Viscosity of Rigid Water Models. J. Chem. Phys. 2010, 132, 096101. (101) Fanourgakis, G. S.; Medina, J.; Prosmiti, R. Determining the Bulk Viscosity of Rigid Water Models. J. Phys. Chem. A 2012, 116, 2564−2570. (102) Song, Y.; Dai, L. L. The Shear Viscosities of Common Water Models by Non-Equilibrium Molecular Dynamics Simulations. Mol. Simul. 2010, 36, 560−567. (103) Zubillaga, R. A.; Labastida, A.; Cruz, B.; Martínez, J. C.; Sánchez, E.; Alejandre, J. Surface Tension of Organic Liquids Using the OPLS/AA Force Field. J. Chem. Theory Comput. 2013, 9, 1611− 1615. (104) Fischer, N. M.; van Maaren, P. J.; Ditz, J. C.; Yildirim, A.; van der Spoel, D. Properties of Organic Liquids When Simulated with Long-Range Lennard-Jones Interactions. J. Chem. Theory Comput. 2015, 11, 2938−2944. (105) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. Dipole Moment of Water from Stark Measurements of H2O, HDO, and D2O. J. Chem. Phys. 1973, 59, 2254−2259. (106) van der Spoel, D.; van Maaren, P. J. The Origin of Layer Structure Artifacts in Simulations of Liquid Water. J. Chem. Theory Comput. 2006, 2, 1−11. (107) Paluch, M. Surface Potential at the Water−Air Interface. Ann. Univ. Mariae Curie-Sklodowska, Sect. AA: Chem. 2015, 70 (2), 1−11. (108) Brodskaya, E.; Zakharov, V. Computer Simulation Study of the Surface Polarization of Pure Polar Liquids. J. Chem. Phys. 1995, 102, 4595−4599. (109) Wilson, M. A.; Pohorille, A.; Pratt, L. R. Comment on ‘‘Study on the Liquid−Vapor Interface of Water. I. Simulation Results of Thermodynamic Properties and Orientational Structure’’. J. Chem. Phys. 1989, 90, 5211−5213. (110) Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B. Affinities of Amino Acid Side Chains for Solvent Water. Biochemistry 1981, 20, 849−855. (111) Sitkoff, D.; Sharp, K. A.; Honig, B. Accurate Calculation of Hydration Free Energies Using Macroscopic Solvent Models. J. Phys. Chem. 1994, 98, 1978−1988. (112) Chakrabarti, P. Geometry of Interaction of Metal Ions with Histidine Residues in Protein Structures. Protein Eng., Des. Sel. 1990, 4, 57−63. (113) Jensen, K. P. Improved Interaction Potentials for Charged Residues in Proteins. J. Phys. Chem. B 2008, 112, 1820−1827. (114) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (115) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The Resp Model. J. Phys. Chem. 1993, 97, 10269−10280. (116) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (117) Cerutti, D. S.; Rice, J. E.; Swope, W. C.; Case, D. A. Derivation of Fixed Partial Charges for Amino Acids Accommodating a Specific Water Model and Implicit Polarization. J. Phys. Chem. B 2013, 117, 2328−2338.

(118) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (119) Jiang, Y.; Zhang, H.; Tan, T. Rational Design of MethodologyIndependent Metal Parameters Using a Nonbonded Dummy Model. J. Chem. Theory Comput. 2016, 12, 3250−3260. (120) Li, P.; Merz, K. M., Jr Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10, 289−297. (121) Duarte, F.; Bauer, P.; Barrozo, A.; Amrein, B. A.; Purg, M.; Åqvist, J.; Kamerlin, S. C. L. Force Field Independent Metal Parameters Using a Nonbonded Dummy Model. J. Phys. Chem. B 2014, 118, 4351−4362. (122) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; HeadGordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549−2564. (123) Lopes, P. E.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D., Jr Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430−5449. (124) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D., Jr An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 116, 4983−5013.

1052


Force Field Benchmark of Amino Acids: I ... - ACS Publications

Recommend Documents