Letter pubs.acs.org/JPCL
Microscopic Origin of Hysteresis in Water Sorption on Protein Matrices Sang Beom Kim, Evan M. Sparano, Rakesh S. Singh, and Pablo G. Debenedetti* Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States S Supporting Information *
ABSTRACT: Despite the importance of water sorption isotherms for a fundamental understanding of protein−water interactions, the microscopic origin of hysteresis between the adsorption and desorption branches is not well understood. Using our recently developed simulation technique, we compute the water sorption isotherms of two proteins, lysozyme and Trp-cage, a miniprotein. We explicitly compare protein−water interactions in adsorption and desorption processes, by analyzing local hydration in terms of hydrogen bonding, water density, and solvent-accessible surface area. We find that significant differences in hydration behavior between adsorption and desorption manifest themselves at the individual amino acid level, in particular around polar or charged residues. We confirm this observation by demonstrating that Trp-cage’s hysteresis can be significantly reduced by mutating charged residues to alanine, a neutral and nonpolar amino acid.
I
the underlying protein−water interactions is still missing. In particular, no clear explanation for the cause of the observed hysteresis has been provided,9,10,14 although several interpretations have been proposed, such as capillary condensation and conformational changes between adsorption and desorption that correspond to the protein being trapped in structurally different basins of its complex underlying energy landscape.1,17 A challenge in understanding protein−water interactions in the context of sorption isotherms has been the lack of appropriate spatiotemporal resolution in experiments, and the absence of efficient simulation techniques to compute the water sorption isotherms of flexible protein matrices. Recent advances in computer simulation methods,9,10 however, have enabled the calculation of sorption isotherms and the consequent systematic investigation of hydration-dependent properties of flexible protein matrices. In this work, we have used the simulation technique of Kim et al.10 to study the microscopic origin of hysteresis in water sorption isotherms through a detailed investigation of local structural and energetic changes along the adsorption and desorption branches of the isotherm. The water sorption isotherms of lysozyme and Trp-cage miniprotein powders were computed at 300 K and 1 bar. The structures of the proteins were taken from the RCSB Protein Data Bank (PDB IDs: 1AKI (lysozyme)18 and 1L2Y (Trpcage)19). Lysozyme has been widely studied on account of its important role in antibacterial resistance and the fact that its size is adequate to be investigated using both the experiments and simulations,20,21 while the Trp-cage miniprotein has been used in numerous computational studies because of its small
n order for proteins to function properly, they must display appropriate structure and dynamics, both of which are profoundly influenced by the interactions between a protein and its hydration water.1−3 Accordingly, proteins must be adequately hydrated under physiological conditions. In low moisture conditions, protein motions are hindered and native structures are disrupted.4−6 For example, the enzymatic activity of lysozyme is significantly reduced as the molecule is dried to below a hydration level of 0.2 g of water/g of dry protein (g/ g).1,7 It is therefore important to understand the interaction between a protein and its surrounding hydration water, especially at hydration levels below approximately 0.4 g/g, above which hydration water begins exhibiting the characteristics of bulk water.1 One of the most widely employed approaches for obtaining insights into protein−water interactions is to study water sorption isotherms on protein powders. The water sorption isotherm provides the relationship between relative humidity at a given temperature and pressure and equilibrium moisture content (typically ranging between 0 to 0.4 g/g). The water sorption isotherms of globular proteins at ambient conditions show pervasive commonalities: sigmoidal shape (type II in the IUPAC convention8) with pronounced hysteresis between adsorption and desorption isotherms.3,9,10 Investigation of water sorption isotherms can be useful in the development of strategies for preservation of protein-based products (prediction of changes to protein stability based on moisture level11), and for lyophilization (prevention of damages to labile biomolecules due to the exposure to severe dehydration during this process12,13). Despite the fact that there have been numerous experimental and computational studies of water sorption isotherms on proteins,3,9,10,14−16 fundamental microscopic understanding of © XXXX American Chemical Society
Received: January 24, 2017 Accepted: February 24, 2017 Published: February 24, 2017 1185
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
Letter
The Journal of Physical Chemistry Letters size, fast folding kinetics, and the fact that it possesses nontrivial secondary and tertiary structures (e.g., α helix and hydrophobic caging).19,22−25 The size range of these proteins (lysozyme, 129 residues; Trp-cage, 20 residues) allows us to probe one aspect of the generality of our findings. The CHARMM2726 and SPC/ E27 force fields were used to model protein and water molecules, respectively, with the Lorentz−Berthelot combining rules. The preparation of the powder systems and simulation of adsorption and desorption were performed by following the procedures described in Kim et al.10 (see Supporting Information for details). In brief, multiple protein units (8 for lysozyme and 16 for Trp-cage) were placed in a simulation box with random translation and rotation. The system was then solvated and equilibrated for 500 ns with NPT molecular dynamics (MD) simulation. The protein powder was dehydrated to h = 0.04 g/g, which is roughly the extent of hydration in freeze-dried protein powders due to residual water.13 Then the protein matrices were hydrated up to 0.4 g/g (adsorption), above which the hydration water begins exhibiting the characteristics of bulk water, and then dehydrated back to 0.04 g/g (desorption). Configurations were saved every h = 0.02 g/g during the adsorption and desorption processes, and each of them was relaxed for 700 ns with NPT MD simulation, of which the last 600 ns was used for data analysis. The water sorption isotherms show the hydration level as a function of relative humidity (RH). RH is defined as RH(%) = 100 × exp((μw − μsat w )/(RT)), where μw is the chemical potential of water, R is the gas constant, and T is temperature. The chemical potential of saturated water vapor (μsat w ) at 300 K was calculated using the ideal gas law and the vapor pressure of SPC/E water (0.01 bar).28 μw was computed using Bennett’s acceptance ratio method29 with numerous test water-molecule insertions and deletions, where the number of insertions and deletions were 20 × 106 and the number of water molecules in the system, respectively.10 The computed water sorption isotherms of lysozyme and Trp-cage (Figure 1) exhibit type II-like shapes (IUPAC convention8) with a pronounced hysteresis between adsorption and desorption isotherms, which is consistent with the typical shape of isotherms of globular proteins. The simulated isotherm of lysozyme shows good agreement with the experimental measurements at ambient temperature16,30−32 (see Figure S1 in the Supporting Information). The difference in RH between adsorption and desorption processes at the same hydration level is largest at approximately h = 0.22 g/g. We therefore compared configurations along the adsorption and desorption branches at this hydration level, in order to identify the differences in protein−water interaction that contribute to the observed hysteresis. We used three metrics to characterize the degree of local hydration: the number of hydrogen bonds (H-bonds) formed between water and individual residues, the number of water molecules within the first hydration shell of the each residue, and the solvent-accessible surface area (SASA) of each residue. A hydrogen bond was considered formed when the donor− acceptor distance is less than 0.35 nm and the acceptor−donorhydrogen angle is less than 30 degrees, where we used the nitrogen and oxygen atoms as acceptors and OH and NH groups as donors.33 A water molecule was considered to be in the first hydration shell of the protein if its oxygen atom is located within a cutoff distance from any non-hydrogen protein atom. This cutoff distance was determined to be 0.31 nm from the radial distribution function between protein non-hydrogen
Figure 1. Water sorption isotherms for lysozyme (top) and Trp-cage (bottom) powders. The orange dashed line indicates the hydration level of 0.22 g/g, at which explicit comparisons between adsorption and desorption were made since at this hydration level the difference in RH is the largest. The black lines are drawn as a guide to the eye.
atoms and water oxygen atoms (see Figure S2 in the Supporting Information). The previous study of Kim et al.10 showed that the adsorption and desorption processes do not exhibit significant differences in these metrics when analyzed for the whole protein. However, this absence of a global signature does not necessarily eliminate the possibility of localized changes in the hydration behavior. Accordingly, we performed the present analysis at a local, microscopic level. We first analyzed the degree of local hydration for each type of amino acid (Figure 2). The number of protein−water H-bonds, number of water molecules within the first hydration shell of each residue, and SASA were first computed for each amino acid, and averaged over the multiple protein units in the computational cell (8 for lysozyme and 16 for Trp-cage). Each quantity was further averaged if more than a single residue is present for each type of amino acid. In general, the differences are rather small and insignificant, suggesting that water molecules hydrate each amino-acid type equally on average between adsorption and desorption processes. However, some amino-acid types (e.g., Glu (E) in lysozyme and Lys (K) in Trp-cage) exhibit significant differences, and these residues have a common trait of having a net charge. This suggests that the chemical nature of amino acids can potentially contribute to the differences in their local hydration between adsorption and desorption processes. We then performed the analysis on each residue basis. As described above, we first computed the degree of local hydration for each residue in the protein matrix, but for this subsequent analysis we averaged only over the multiple 1186
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
Letter
The Journal of Physical Chemistry Letters
Figure 2. Comparison between the adsorption and desorption in terms of (A) the number of protein−water H-bonds, (B) number of water molecules within the first hydration shell, and (C) solvent-accessible surface area (SASA) for each type of amino acid in the lysozyme and Trp-cage powders.
exhibits a large SASA difference between adsorption and desorption. This figure clearly shows that the ASP9 residue is significantly more buried by the nearby residues (shown in gray surface representation) during the desorption process, supporting the importance of localized differences in protein hydration. In order to further understand which type of amino acid contributes more to these localized measures of hydration hysteresis, in Figure 3C we plotted the magnitude (the absolute value) of difference in the number of protein−water H-bonds between adsorption and desorption as a function of amino-acid type (averaged if more than a single residue is present for each amino-acid type). When the amino acid types are arranged in ascending order of the differences, it can be seen that the largest differences occur predominantly for charged and polar amino acids. This contrast between the behavior of apolar and polar or charged residues is more significant for lysozyme than for Trpcage, perhaps reflecting the former’s larger size and greater chemical heterogeneity. The fact that the charged residues exhibit more significant differences in local hydration between adsorption and desorption than neutral ones may suggest that the former are largely responsible for the hysteresis. In order to confirm this hypothesis, each of Trp-cage’s charged residues (K8, D9, and
proteins, not over multiple occurrences of the same amino acid in the protein. For example, 9ALA and 11ALA were averaged together in order to compute the ALA-specific contributions in Figure 2, but they are now considered separately. Figure 3A shows the analysis for lysozyme and Trp-cage, where each residue was colored by the difference between the adsorption and desorption values of the degree of local hydration, measured by the number of protein−water H-bonds, the number of water molecules in the first hydration shell, and SASA. The red or blue colors indicate that these residues are more hydrated during adsorption or desorption, respectively. This suggests that water molecules interact preferentially with different sites when they are being adsorbed on, or desorbed from, the protein. This can potentially cause a difference in the chemical potential of water (hence the relative humidity) upon adsorption versus upon desorption, even at the same hydration level. Moreover, this difference is more pronounced for lysozyme than Trp-cage, consistently with the larger protein size, and the consequently greater heterogeneity. Figure 3B shows representative configurations of Trp-cage from the adsorption and desorption processes, following structural alignment using the backbone atoms.34,35 The side chain of the ASP9 residue is shown explicitly because this residue 1187
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
Letter
The Journal of Physical Chemistry Letters
Figure 3. (A) Graphical representation of the local hydration of lysozyme (top) and Trp-cage (bottom). A residue was colored red if it is more locally hydrated during adsorption, and blue if it is more locally hydrated during desorption. White residues exhibit no difference between the adsorption and desorption processes. (B) Representative configurations of Trp-cage during adsorption and desorption processes. Aspartic acid (9ASP) is explicitly shown, while the rest of the protein is represented with the backbone and gray solvent-accessible surface. (C) The absolute value of the difference in the number of protein−water H-bonds between adsorption and desorption as a function of amino acid type. The largest differences occur mostly for charged (orange) and polar (green) residues.
elucidate the generality of our results based so far on two proteins.
R16) was mutated to alanine, a nonpolar amino acid. Trp-cage was chosen for this mutation study because of the comparatively smaller computational cost associated with its smaller size. The water sorption isotherm of the mutated Trpcage (Figure 4A) indeed shows significantly reduced adsorption−desorption hysteresis. When the local hydration level is visualized (Figure 4B), we see a considerable reduction in the local differences in hydration between adsorption and desorption when compared to the unmutated Trp-cage (Figure 3A). In summary, we have identified the difference in local indicators of hydration between adsorption and desorption processes as the microscopic signature of hysteresis in the water sorption isotherm for lysozyme and Trp-cage powders. The local hydration differences occur at an individual residue level, and are more strongly localized at charged and polar residues. This was confirmed by mutating the charged amino acids in Trp-cage to alanine, which significantly reduced the magnitude of the hysteresis. Our results suggest that localized differences in hydration can cause correspondingly significant differences in the relative humidity that the protein matrix equilibrates to at a given hydration level, despite the fact that there may not be notable differences in hydration at the whole-protein level, or after averaging over different occurrences of amino-acid types. A potentially interesting avenue of future inquiry will be to investigate the local hydration behavior of structurally different proteins using a similar computational approach, in order to
■
COMPUTATIONAL METHODS
For all MD simulations, we used the GROMACS36,37 simulation package. Equations of motion were integrated using the leapfrog algorithm with a 2 fs time step. The Nosé−Hoover thermostat38,39 with a 0.2 ps time constant and the Parrinello−Rahman barostat40,41 with a 2 ps time constant were respectively used to maintain the temperature and pressures at 300 K and 1 bar. Anisotropic pressure coupling was implemented in order to allow independent fluctuations of the dimensions of the simulation box. Periodic boundary conditions were applied in all three dimensions. The shortrange interactions were truncated at 1 nm, and the standard long-range dispersion corrections were applied for the pressure and energy.42 The smooth particle mesh Ewald method43 was used to compute the reciprocal part of the Ewald sum for the long-range electrostatics. All bonds in the protein and water molecules were constrained using the linear constraint solver algorithm (LINCS)44,45 and SETTLE,46 respectively. Each simulation trajectory was divided into five equally sized blocks (120 ns) in order to compute the standard errors using the block-averaging analysis.47 1188
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
Letter
The Journal of Physical Chemistry Letters
(2) Teeter, M. M. Water-Protein Interactions: Theory and Experiment. Annu. Rev. Biophys. Biophys. Chem. 1991, 20, 577−600. (3) Panagopoulou, A.; Kyritsis, A.; Aravantinou, A.-M.; Nanopoulos, D.; i Serra, R. S.; Ribelles, J. L. G.; Shinyashiki, N.; Pissis, P. Glass Transition and Dynamics in Lysozyme-Water Mixtures Over Wide Ranges of Composition. Food Biophys. 2011, 6, 199−209. (4) Prestrelski, S. J.; Tedeschi, N.; Arakawa, T.; Carpenter, J. F. Dehydration-Induced Conformational Transitions in Proteins and Their Inhibition by Stabilizers. Biophys. J. 1993, 65, 661−671. (5) Roh, J. H.; Curtis, J. E.; Azzam, S.; Novikov, V. N.; Peral, I.; Chowdhuri, Z.; Gregory, R. B.; Sokolov, A. P. Influence of Hydration on the Dynamics of Lysozyme. Biophys. J. 2006, 91, 2573−2588. (6) Kim, S. B.; Gupta, D. R.; Debenedetti, P. G. Computational Investigation of Dynamical Transitions in Trp-Cage Miniprotein Powders. Sci. Rep. 2016, 6, 25612. (7) Lind, P. A.; Daniel, R. M.; Monk, C.; Dunn, R. V. Esterase Catalysis of Substrate Vapour: Enzyme Activity Occurs at Very Low Hydration. Biochim. Biophys. Acta, Proteins Proteomics 2004, 1702, 103−110. (8) Sing, K. S. W.; Everett, D. H.; Haul, R. A. W.; Moscou, L.; Pierotti, R. A.; Rouquerol, J.; Siemieniewska, T. Reporting Physisorption Data for Gas Solid Systems with Special Reference to the Determination of Surface-Area and Porosity (Recommendations 1984). Pure Appl. Chem. 1985, 57, 603−619. (9) Palmer, J. C.; Debenedetti, P. G. Computer Simulation of Water Sorption on Flexible Protein Crystals. J. Phys. Chem. Lett. 2012, 3, 2713−2718. (10) Kim, S. B.; Palmer, J. C.; Debenedetti, P. G. A Computational Study of the Effect of Matrix Structural Order on Water Sorption by Trp-Cage Miniproteins. J. Phys. Chem. B 2015, 119, 1847−1856. (11) Bell, L. N.; Labuza, T. P. Moisture Sorption: Practical Aspects of Isotherm Measurement and Use, 2nd ed; American Association of Cereal Chemists: St. Paul, MN, 2000. (12) Carpenter, J. F.; Pikal, M. J.; Chang, B. S.; Randolph, T. W. Rational Design of Stable Lyophilized Protein Formulations: Some Practical Advice. Pharm. Res. 1997, 14, 969−975. (13) Wang, W. Lyophilization and Development of Solid Protein Pharmaceuticals. Int. J. Pharm. 2000, 203, 1−60. (14) Teng, C. D.; Zarrintan, M. H.; Groves, M. J. Water Vapor Adsorption and Desorption Isotherms of Biologically Active Proteins. Pharm. Res. 1991, 8, 191−195. (15) Kumagai, H.; Seto, H.; Sakurai, H.; Ishii, K.; Kumagai, H. Analysis of Water Sorption Behavior of Native and Denatured Proteins by Solution Thermodynamics. Biosci., Biotechnol., Biochem. 1997, 61, 1307−1311. (16) Lű scher-Mattli, M.; Rű egg, M. Thermodynamic Functions of Biopolymer Hydration. I. Their Determination by Vapor Pressure Studies, Discussed in an Analysis of the Primary Hydration Process. Biopolymers 1982, 21, 403−418. (17) Bryan, W. P. Thermodynamic Models for Water−Protein Sorption Hysteresis. Biopolymers 1987, 26, 1705−1716. (18) Artymiuk, P.; Blake, C.; Rice, D.; Wilson, K. The Structures of the Monoclinic and Orthorhombic Forms of Hen Egg-White Lysozyme at 6 Å Resolution. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1982, 38, 778−783. (19) Neidigh, J. W.; Fesinmeyer, R. M.; Andersen, N. H. Designing a 20-Residue Protein. Nat. Struct. Biol. 2002, 9, 425−430. (20) Mckenzie, H. A.; White, F. H. Lysozyme and α-Lactalbumin: Structure, Function, and Interrelationships. Adv. Protein Chem. 1991, 41, 173−315. (21) Davies, D. R.; Cohen, G. H. Interactions of Protein Antigens with Antibodies. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 7−12. (22) Kim, S. B.; Palmer, J. C.; Debenedetti, P. G. Computational Investigation of Cold Denaturation in the Trp-cage Miniprotein. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 8991−8996. (23) Kim, S. B.; Dsilva, C. J.; Kevrekidis, I. G.; Debenedetti, P. G. Systematic Characterization of Protein Folding Pathways using Diffusion Maps: Application to Trp-Cage Miniprotein. J. Chem. Phys. 2015, 142, 085101.
Figure 4. (A) Water sorption isotherm for mutated Trp-cage powder. (B) Graphical representation of the local hydration of mutated Trpcage, in terms of the number of water molecules in the first hydration shell and SASA. The color scheme is the same as Figure 3A.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b00184. (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Pablo G. Debenedetti: 0000-0003-1881-1728 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS P.G.D. gratefully acknowledges financial support from the National Science Foundation (Grant No. CBET-1263565). The computations were performed at the Terascale Infrastructure for Groundbreaking Research in Engineering and Science (TIGRESS), at Princeton University. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (Grant No. ACI-1053575).
■
REFERENCES
(1) Rupley, J. A.; Careri, G. Protein Hydration and Function. Adv. Protein Chem. 1991, 41, 37−172. 1189
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
Letter
The Journal of Physical Chemistry Letters (24) Paschek, D.; Hempel, S.; García, A. E. Computing the Stability Diagram of the Trp-Cage Miniprotein. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 17754−17759. (25) Juraszek, J.; Bolhuis, P. G. Sampling the Multiple Folding Mechanisms of Trp-Cage in Explicit Solvent. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 15859−15864. (26) MacKerell, A. D.; Feig, M.; Brooks, C. L. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400−1415. (27) Berendsen, H.; Grigera, J.; Straatsma, T. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (28) Vega, C.; Abascal, J. L. F.; Nezbeda, I. Vapor-Liquid Equilibria from the Triple Point up to the Critical Point for the New Generation of TIP4P-Like Models: TIP4P/Ew, TIP4P/2005, and TIP4P/ice. J. Chem. Phys. 2006, 125, 034503. (29) Bennett, C. H. Efficient Estimation of Free-Energy Differences from Monte-Carlo Data. J. Comput. Phys. 1976, 22, 245−268. (30) Hnojewyj, W. S.; Reyerson, L. H. Further Studies on the Sorption of H2O and D2O Vapors by Lysozyme and the DeuteriumHydrogen Exchange Effect. J. Phys. Chem. 1961, 65, 1694−1698. (31) Smith, A. L.; Shirazi, H. M.; Mulligan, S. R. Water Sorption Isotherms and Enthalpies of Water Sorption by Lysozyme Using the Quartz Crystal Microbalance/Heat Conduction Calorimeter. Biochim. Biophys. Acta, Protein Struct. Mol. Enzymol. 2002, 1594, 150−159. (32) D’Arcy, R. L.; Watt, I. C. Analysis of Sorption Isotherms of Non-Homogeneous Sorbents. Trans. Faraday Soc. 1970, 66, 1236− 1245. (33) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, N. Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media. J. Phys. Chem. B 2006, 110, 4393−4398. (34) Kabsch, W. A Discussion of the Solution for the Best Rotation to Relate Two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1978, 34, 827−828. (35) Kabsch, W. A Solution for the Best Rotation to Relate Two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1976, 32, 922−923. (36) 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. (37) Berendsen, H. J.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (38) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (39) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (40) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular-Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (41) Nosé, S.; Klein, M. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (42) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (43) 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. (44) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (45) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (46) 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. (47) Flyvbjerg, H.; Petersen, H. G. Error Estimates on Averages of Correlated Data. J. Chem. Phys. 1989, 91, 461−466.
1190
DOI: 10.1021/acs.jpclett.7b00184 J. Phys. Chem. Lett. 2017, 8, 1185−1190
