Enhanced Conformational Sampling of Peptides via Reduced Side

Nov 15, 2010 - ... the spirit of Bennett's mass-tensor dynamics ( Bennett , C. H. J. Comput. ... Chad W. Hopkins , Scott Le Grand , Ross C. Walker , a...
0 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 15935–15940

15935

Enhanced Conformational Sampling of Peptides via Reduced Side-Chain and Solvent Masses I-Chun Lin Department of Chemistry, New York UniVersity, New York, New York 10003, United States

Mark E. Tuckerman* Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York 10003, United States ReceiVed: October 14, 2010

When only equilibrium thermodynamics or averages of position-dependent functions in classical molecular dynamics calculations are sought, the choice of the atomic masses in the Hamiltonian becomes irrelevant. Consequently, the masses can be used as free parameters to help improve conformational sampling efficiency in the spirit of Bennett’s mass-tensor dynamics (Bennett, C. H. J. Comput. Phys. 1975, 19, 267.). Here we report on the effect of reducing side-chain and solvent masses on the folding behavior of a 9-residue hairpinforming peptide in solution. A physical motivation for reducing the solvent mass is an effective reduction in solvent viscosity. Side-chain mass scaling is motivated by the idea of solvent potentials of mean force which are employed in select coarse-grained protein models. In the limit of very large mass differences, the mass reduction effectively creates an adiabatic decoupling between solvent, side-chain, and backbone motions, so that both the backbone and side chains move on the instantaneous solvent potential of mean force (PMF) surface, and the backbone additionally moves on the side-chain PMF. Because of the arbitrariness in the choice of masses, this limit only needs to be reached approximately in practice. In particular, we show that a 10-fold reduction in solvent masses and a side-chain mass scale that is intermediate between the scaled solvent and the backbone lead to a quantitative enhancement in conformational sampling. 1 H ) pTM-1p + V(r) 2

1. Introduction A correct three-dimensional structure is essential for a protein to function properly,1 and various techniques have been proposed to predict a protein’s tertiary structure from its amino acid sequence alone. Nevertheless, sampling of the statistical ensemble distribution based on current computational methods is often insufficient to probe the conformational space of systems described by complex, multidimensional energy landscapes that are typical of proteins. This difficulty has limited the impact of computer simulations on molecular biology, in particular protein structure prediction and design. Computational studies of protein structures fall within a class of calculations aimed at predicting the equilibrium properties of a system, and these are obtained from averages of phase-space functions that depend only on the atomic positions. When such averages are sought in a canonical (NVT) or an isothermal-isobaric (NPT) simulation, the choice of the individual atomic masses in the Hamiltonian becomes irrelevant. These masses can thus be used as free parameters to help improve the efficiency of conformational sampling. Moreover, as was recognized some time ago by Bennett,2 the usual diagonal mass matrix could be replaced by a full mass tensor whose elements could be chosen so as to bring high- and low-frequency motions onto a similar time scale. If p and r represent the full set of momenta and coordinates in a system, then Bennett’s mass-tensor dynamics is tantamount to employing a Hamiltonian of the form * To whom correspondence should be addressed. E-mail: mark.tuckerman@ nyu.edu. Tel: +1 212 9988471. Fax: +1 212 2607905.

(1)

where V(r) is the potential, and M-1 is the inverse of a generally nondiagonal mass tensor. Clearly, even for the Hamiltonian in eq 1, the canonical average of a position-dependent function O(r)

∫ dp e-βp M p/2 ∫ dr O(r)e-βV(r) ∫ dp e-βp M p/2 ∫ dr e-βV(r) ∫ dr O(r)e-βV(r) ∫ dr e-βV(r) T

〈O(r)〉 )

)

-1

T

-1

(2)

is independent of the choice of the mass tensor. The idea of choosing masses in order to reduce the range of time scales has proved to be highly useful in path integral molecular dynamics,3 for example. Others have used mass scaling to damp high angular velocities of water molecules in simulations4 and to increase the effective momentum of dihedral rotational motion in barrier crossing.5 As a concrete example, consider reducing the mass of the solvent by an amount sufficient to create an adiabatic decoupling between the solvent and the solute. In a spirit similar to adiabatic free-energy methods,6-9 the mass disparity would cause the solute to move on the instantaneous solvent potential of mean force (PMF) surface. In fact, weighting the masses of solvent molecules in explicit-solvent simulations by λ changes the

10.1021/jp109865y  2010 American Chemical Society Published on Web 11/15/2010

15936

J. Phys. Chem. B, Vol. 114, No. 48, 2010

Lin and Tuckerman

TABLE 1: Average Total and Potential Energies (Etot and EPE, in kcal/mol), Temperatures (T, in K) and Pressures (P, in atm), and the Energy Drift Per Degree of Freedom (in kBT/ns) during Run A45 a Etot EPE T P Edrift a

control

solv0.1

solv0.1-side0.8

solv0.1-side0.6

solv0.1-side0.4

-22382.6 (91.6) -27402.1 (73.4) 300.0 (3.3) -26 (247) 0.6 × 10-4

-22382.0 (91.6) -27400.8 (73.4) 300.0 (3.3) -27 (248) 1.6 × 10-4

-22382.6 (91.6) -27401.4 (72.8) 300.0 (3.3) -26 (248) 2.7 × 10-4

-22382.0 (92.9) -27401.4 (74.0) 300.0 (3.3) -27 (247) 1.8 × 10-4

-22382.6 (91.7) -27401.9 (73.6) 300.0 (3.3) -27 (247) 1.1 × 10-4

Standard deviations are given in parentheses.

effective viscosity by λ.10,11 It is well-known that solvent viscosity influences the dynamics of protein folding: both experimental and theoretical molecular dynamics (MD) studies have shown that proteins exhibit more flexibility and shorter correlation times in low-viscosity solvents.12-21 A decrease in viscosity thus allows the folding to accelerate until other phenomena (e.g., intrachain interactions) take over, and by artificially lowering the effective solvent viscosity in simulations, more efficient sampling can be achieved.14,16,20 Obviously, effective viscosity is irrelevant for thermodynamic properties as long as the fundamental solvent-solvent and solvent-solute interactions remain unaltered and the Boltzmann distribution is sufficiently sampled.14 Nevertheless, this discussion gives a physical motivation for using mass scaling to study static equilibrium properties of solvated peptides. The idea of mass scaling can be applied beyond just reducing the time scale of solvent motion. In particular, the motion of side chains also influences sampling efficiency. For this reason, some coarse-grained force fields, such as the latest version of UNRES by Scheraga and co-workers,22 employ rotamer PMFs as part of the coarse-grained potential energy function. Here, we propose that reducing the masses of side-chain atoms, together with a reduction in solvent masses, can be highly effective in improving conformational sampling for all-atom force-field models. If the side chains are adiabatically decoupled from the backbone, then the backbone will move on the instantaneous side-chain or rotamer PMF surface for the force field employed. We argue further that the time scale of sidechain motion should be intermediate between that of solvent and backbone so that the side chains remain adiabatically decoupled from the solvent and, therefore, move on the instantaneous solvent PMF surface. In practice, sampling enhancement can be achieved without fully reaching the ideal adiabatically decoupled limit. To demonstrate the utility of the proposed hierarchy of mass scaling, we sampled the conformational space of a solvated peptide Tyr-Gln-Asn-Pro-Asp-Gly-Ser-Gln-Ala (YQNPDGSQA)23 at 300 K under an ambient pressure of 1 atm. When folded, this mutated peptide has enhanced turn-forming propensities and exhibits NOEs compatible with a β-hairpin structure.24,25 Its small size makes simulations with explicit water not too expensive to carry out and has been the subject of several computational studies.26-30

points for subsequent runs. In the “control” simulation, solvent and side-chain masses were unaltered. In the “solv0.1” simulation, all atoms of TIP3P waters had masses that were one-tenth of the original to mimic a low-viscosity solvent (the new viscosity is approximately 0.3 times the original). Finally, in a third set of simulations, water masses were scaled by a factor of 0.1, and the masses of side-chain atoms in all amino acids (except proline and glycine) were scaled by 0.8, 0.6, and 0.4, corresponding to simulations “solv0.1-side0.8”, “solv0.1-side0.6”, and “solv0.1-side0.4”, respectively. For better statistics, three independent runs (runs A, B, and C) starting from different initial configurations and velocities were carried out for each simulation setup. The fully extended peptide was set up using PepBuild.32 All simulations were carried out with the PINY_MD code33 and the CHARMM22 force field.34,35 Periodic boundary conditions were applied throughout. The SHAKE36 and RATTLE37 algo-

2. Computational Details In order to prepare the system, a fully extended 9-residue peptide and one sodium counterion were solvated in a cubic box of 2743 TIP3P31 waters. The system was equilibrated in the canonical (NVT) ensemble for 100 ps, then in the isothermal-isobaric (NPT) ensemble for 100 ps, leading to an average box length of 43.5 Å, and again in the NVT ensemble for another 100 ps. Configurations from the equilibrated system (see Figure 4a for one example) were then used as the starting

Figure 1. Two-dimensional probability distributions as a function of the radius of gyration (Rgyr) and the CR rmsd from the native state. Data are collected over three independent runs, each of 2.8 × 107 simulation steps.

Enhanced Conformational Sampling of Peptides

J. Phys. Chem. B, Vol. 114, No. 48, 2010 15937

Figure 2. Time evolution of the CR (black) and heavy atom (red) rmsds from the native β-hairpin for the (a) control, (b) solv0.1, (c) solv0.1-side0.8, (d) solv0.1-side0.6, and (e) solv0.1-side0.4 setups. A red dashed line is drawn at 1 Å.

rithms were used to constrain the water OH bonds to relative tolerances of 10-6 and 10-5 for the NVT and NPT runs, respectively. The electrostatics were treated with the Ewald summation technique,38 and the short-range forces were switched off at 10 Å. The Nose´-Hoover chain (NHC) method39 was used for the thermostating (coupling time 160 fs, chain length of 2, Suzuki/Yoshida order of 3, and decomposition number of 2). Each degree of freedom was coupled to an individual thermostat. The coupling time for the barostat was 600 fs. The equations of motion were integrated using the multiple time scale algorithm, r-RESPA,40,41 to exploit the various separations of time scales. Two r-RESPA levels were used: long-range forces were integrated using two steps per full MD step, and all intramolecular forces were integrated using eight steps per longrange force step. The outermost time step was 0.25 fs for all simulations involved mass scaling and 0.50 fs for the rest. The time step was halved to account for the higher frequency induced by the mass scaling. The small time step here is to ensure a tight energy conservation (see Table 1) and can be further relaxed, especially if bonds involving hydrogen are constrained. All simulations were carried out for 2.8 × 107 steps, corresponding to 14 ns for the control and 7 ns for the rest, and atomic coordinates were saved every 500 steps. Figures 3 and 4 were rendered with VMD.42

Figure 3. Structures with small CR rmsds from the solv0.1-side0.8 setup in (a) run A, (b) run B, and (c) run C. The corresponding CR rmsd and radius of gyration (in Å) are listed below.

3. Results and Discussion As mentioned previously, the choice of particle masses does not affect any thermodynamic property or equilibrium average of a position-dependent function. Some examples of quantities whose averages do not depend on the choice of masses (average energies, temperatures, and pressures) are tabulated in Table 1, and the differences in these quantities over the various simulations are negligible. We used the radius of gyration (Rgyr) of the heavy backbone atoms (excluding oxygen) and the root-mean-square deviation (rmsd) from the native β-hairpin (a folded structure from ref 30 is used as the reference state) to track conformational changes that the peptide undergoes. The rmsd describes how similar the instantaneous structure is to the native state, while Rgyr gives a measure of the peptide’s overall size. The peptide is considered folded if its CR rmsd is less than 1 Å from the native state reported in ref 30 (comparisons to the one of ref 28 give similar results). In a real biological system, proteins exist in an equilibrium between unfolded and folded states. Even though our simulations are not long enough to observe this continuous

Figure 4. (a) Initial and a (b) folded structures of Tyr-Gln-Asn-ProAsp-Gly-Ser-Gln-Ala from the solv0.1-side0.6 setup in run A. Interstrand H-bonds are depicted in (b).

reversible folding process in all setups, Figures 1 and 2 demonstrate the improvement in sampling by mass scaling. With mass scaling, the peptide visits more of the available configuration space. Small rmsd and Rgyr are obtained when solvent masses only are scaled, but none manages to reach the folded state. One implicit-solvent MD study16 has shown that for solvent viscosity down to just above a tenth of the water viscosity, the folding rate depends linearly on the viscosity and

15938

J. Phys. Chem. B, Vol. 114, No. 48, 2010

Lin and Tuckerman

Figure 5. Values of d1, d2, d3, and d4 during run B. A red dashed line is drawn at 0.25 nm.

can be explained by the reaction rate theory of Kramers.43 Since the viscosity in the solv0.1 setup is only 3 times lower, the length of our simulations is still insufficient for a hairpin of this size to form (the typical folding rate of a β-hairpin is in the order of hundreds of nanoseconds or microseconds44). However, when the side-chain masses are reduced as well, the sampling is improved further. There is a higher probability to visit states with lower CR rmsd (less than 2 Å) in both solv0.1-side0.8 and solv0.1-side0.6 setups. The hairpin backbone structure is (nearly) formed in the solv0.1-side0.8 setup, as the CR rmsds for some structures are close to or even lower than 1 Å (Figure 2). In many cases, however, the corresponding heavy-atom rmsd is much larger, indicating that the positions of side chains still deviate from those of the reference state. Several such structures are depicted in Figure 3. The solv0.1-side0.6 setup, on the other hand, reaches the folded state in both runs A and B, and the peptide stays folded for approximately 2.5 and 1 ns, respectively.

Figure 4b shows one of the folded conformations. Curiously enough, by reducing the side-chain masses even further (solv0.1-side0.4), sampling does not seem to become more efficient. This observation underscores the notion that the mass disparity between side chains and solvent should be large enough to approximate an adiabatic decoupling between side-chain and solvent motions. As a side note, the diffusion coefficient of water in setups where the water mass is scaled down is 3 times as high as the control, in agreement with previous studies.10,11 The folded structure has a turn at positions Asn3-Gly6 and interstrand hydrogen bonds (H-bonds) that are characteristics of a β-hairpin. Three such H-bonds are depicted in Figure 4b: d1 (between O of Asn3 and H of Gly6), d2 (between H of Asn3 and O of Ser7), and d3 (between O of Tyr1 and H of Ala9). Using the same definition as ref 24, a H-bond is considered present if the distance between the donor H and the acceptor O is less than 2.5 Å. An additional distance, d4, is defined as the

Enhanced Conformational Sampling of Peptides distance between the Cδ atoms of Gln2 and Gln8. These two residues form the only long-range contact in the native state that is not involved in forming the hairpin turn, and the contact serves as a measure of the time ordering of tertiary and secondary structure formation.30 The time evolution of these distances from run B is shown in Figure 5 as an example.45 The simulation with only reduced water mass shows little improvement, but all setups with reduced side-chain masses manage to visit values that correspond to a typical H-bond length in at least two of the distances. In addition, it is only in the solv0.1-side0.6 setup that we observe simultaneously small values of d1, d2, d3, and d4. Reducing the masses of side chains increases the side-chain motion; this allows the peptide to find its folded structure more readily. Figure 5 also sheds some light on the possible folding mechanism. It was postulated in ref 30 that, at 280 K, d1 is formed first, followed by d2, and finally by d3 (a hairpin zipper model), but more kinetic folding pathways open up at higher temperatures. Our results, however, show a simultaneous decrease in the distances between the interstrand H-bonding atoms, indicating a highly cooperative folding process that is closer to the mechanism proposed by refs 27 and 28. This cooperative behavior is also observed in run A (see Supporting Information). The differences among the studies can be attributed to several factors. First, our simulated temperature is slightly elevated to accelerate the sampling process. Second, the force fields in these studies are all different. Furthermore, the conformational distributions in our simulations are still not close to that characteristic of the “folding-unfolding” equilibrium. Caution thus needs to be taken when interpreting the results, and much longer trajectories would be needed to decipher the underlying mechanism. Nevertheless, the key point demonstrated here is that the mass hierarchy of light solvent, intermediate side-chain, and heavy backbone leads to enhanced conformational sampling efficiency and is a scheme that could be used in conjunction with other enhanced sampling approaches. 4. Conclusion We have shown that reducing the masses of both solvent molecules and side-chain atoms leads to a qualitative enhancement in the efficiency of conformational sampling in an allatom model of a small peptide. Under appropriate conditions of adiabatic decoupling, a slow subsystem can be made to move on the PMF generated by a fast subsystem.46 We have demonstrated here that a hierarchy of solvent, side-chain, and backbone time scales in a peptide can be created using appropriate mass scalings such that the backbone and side chains move on an approximate solvent PMF and the backbone additionally moves on an approximate side-chain PMF. Our results indicate that within the same number of time steps, a simulation with both solvent and side-chain masses scaled is able to visit the folded state of this β-hairpin peptide. It is also observed that choosing side-chain masses too low degrades the improvement in sampling, indicating the importance of having the solvent and side-chain motions approximately adiabatically decoupled. Since the simulation time here is too short to reach conformational equilibria, a more extensive examination on the correlation between the side-chain masses and the improvement in sampling would be needed. One significant challenge that faces MD simulations employing enhanced conformational sampling techniques such as replica-exchange MD47-49 or the reference potential spatial warping algorithm50,51 is the treatment of explicit solvent. Since a peptide will undergo rapid conformational changes under the effect of these techniques, explicit solvent could severely hinder

J. Phys. Chem. B, Vol. 114, No. 48, 2010 15939 rapid conformational changes if it cannot follow the motion of the solute. Mass scaling can thus be used in synergy with these approaches both to further enhance sampling and to ensure that the solvent can follow rapid changes in solute conformations. Acknowledgment. I.-C.L. thanks the Swiss National Science Foundation (PBELP2-123062) for its financial support, and M.E.T. acknowledges support from NSF CHE-0704036 and CHE-1012545. Supporting Information Available: Additional information as noted in text. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Branden, C.; Tooze, J. Introduction to Protein Structure; Garland Publishing: New York, 1999. (2) Bennett, C. H. J. Comput. Phys. 1975, 19, 267. (3) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. J. Chem. Phys. 1993, 99, 2796. (4) Pomes, R.; McCammon, J. A. Chem. Phys. Lett. 1990, 166, 425. (5) Mao, B.; Friedman, A. R. Biophys. J. 1990, 58, 803. (6) Rosso, L.; Minary, P.; Zhu, Z.; Tuckerman, M. E. J. Chem. Phys. 2002, 116, 4389–4402. (7) VandeVondele, J.; Rothlisberger, U. J. Phys. Chem. B 2002, 106, 203–208. (8) Maragliano, L.; Vanden-Eijnden, E. Chem. Phys. Lett. 2006, 426, 168–175. (9) Abrams, J. B.; Tuckerman, M. E. J. Phys. Chem. B 2008, 112, 15742–15757. (10) Walser, R.; Mark, A. E.; van Gunsteren, W. F. Chem. Phys. Lett. 1999, 303, 583–586. (11) Walser, R.; Hess, B.; Mark, A. E.; van Gunsteren, W. F. Chem. Phys. Lett. 2001, 334, 337–342. (12) Jacob, M.; Schindler, T.; Balbach, J.; Schmid, F. X. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 5622–5627. (13) Klimov, D. K.; Thirumalai, D. Phys. ReV. Lett. 1997, 79, 317– 320. (14) Walser, R.; van Gunsteren, W. F. Proteins: Struct., Funct, Genet. 2001, 42, 414–421. (15) Jas, G. S.; Eaton, W. A.; Hofrichter, J. J. Phys. Chem. B 2001, 105, 261–272. (16) Zagrovic, B.; Pande, V. J. Comput. Chem. 2003, 24, 1432–1436. (17) Qiu, L.; Hagen, S. J. J. Am. Chem. Soc. 2004, 126, 3398–3399. (18) Pabit, S. A.; Roder, H.; Hagen, S. J. Biochemistry 2004, 43, 12532– 12538. (19) Hagen, S. J.; Qiu, L.; Pabit, S. A. J. Phys.: Condens. Matter 2005, 17, S1503–S1514. (20) Gee, P. J.; van Gunsteren, W. F. Chem.sEur. J. 2006, 12, 72–75. (21) Rhee, Y. M.; Pande, V. S. J. Phys. Chem. B 2008, 112, 6221– 6227. (22) Kozlowska, U.; Liwo, A.; Scheraga, H. A. J. Comput. Chem. 2010, 31, 1143–1153. (23) Blanco, F. J.; Jimenez, M. A.; Herranz, J.; Rico, M.; Santoro, J.; Nieto, J. L. J. Am. Chem. Soc. 1993, 115, 5887–5888. (24) Constantine, K. L.; Mueller, L.; Andersen, N. H.; Tong, H.; Wandler, C. F.; Friedrichs, M. S.; Bruccoleri, R. E. J. Am. Chem. Soc. 1995, 117, 10841–10854. (25) Friedrichs, M. S.; Stouch, T. R.; Bruccoleri, R. E.; Mueller, L.; Constantine, K. L. J. Am. Chem. Soc. 1995, 117, 10855–10864. (26) Constantine, K. L.; Friedrichs, M. S.; Stouch, T. R. Biopolymers 1996, 39, 591–614. (27) Wu, W.; Wang, S.; Brooks, B. R. J. Am. Chem. Soc. 2002, 124, 5282–5283. (28) Wu, W.; Brooks, B. R. Biophys. J. 2004, 86, 1946–1958. (29) Baumketner, A.; Shea, J.-E. J. Phys. Chem. B 2005, 109, 21322– 21328. (30) Baumketner, A.; Shea, J.-E. Theor. Chem. Acc. 2006, 116, 262– 273. (31) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (32) Singh, B.; Kaur, T. PepBuild, http://www.imtech.res.in/bvs/ pepbuild/. (33) Tuckerman, M. E.; Yarne, D. A.; Samuelson, S. O.; Hughes, A. L.; Martyna, G. J. Comput. Phys. Commun. 2000, 128, 333–376. (34) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. J. Phys. Chem. B 1998, 102, 3586–3616.

15940

J. Phys. Chem. B, Vol. 114, No. 48, 2010

(35) MacKerell, A. D., Jr.,; Feig, M.; Brooks, C. L., III. J. Comput. Chem. 2004, 25, 1400–1415. (36) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (37) Andersen, H. C. J. Comput. Phys. 1983, 52, 24–34. (38) Ewald, P. P. Ann. Phys. 1921, 369, 253–287. (39) Martyna, G. J.; Tuckerman, M. E.; Klein, M. L. J. Chem. Phys. 1992, 97, 2635–2643. (40) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990–2001. (41) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117–1157. (42) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. (43) Kramers, H. A. Physica 1940, 7, 284–304.

Lin and Tuckerman (44) Kubelka, J.; Hofrichter, J.; Eaton, W. A. Curr. Opin. Struct. Biol. 2004, 14, 76. (45) Results from other two runs are included in the Supporting Information. (46) Rosso, L.; Tuckerman, M. E. Mol. Simul. 2002, 28, 91–112. (47) Hukushima, K.; Nemoto, K. J. J. Phys. Soc. Jpn. 1996, 65, 1604– 1608. (48) Hansmann, U. H. E. Chem. Phys. Lett. 1997, 281, 140–150. (49) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141–151. (50) Zhu, Z.; Tuckerman, M. E.; Samuelson, S. O.; Martyna, G. J. Phys. ReV. Lett. 2002, 88, 100201. (51) Minary, P.; Tuckerman, M. E.; Martyna, G. J. SIAM J. Sci. Comput. 2008, 30, 2055–2083.

JP109865Y