ARTICLE pubs.acs.org/JPCB
Enhancing the Configurational Sampling of Ions in Aqueous Solution Using Adiabatic Decoupling with Translational Temperature Scaling Anna-Pitschna E. Kunz and Wilfred F. van Gunsteren* Laboratory of Physical Chemistry, Swiss Federal Institute of Technology ETH, 8093 Z€urich, Switzerland ABSTRACT: Three methods to enhance the configurational sampling of ions in aqueous solution, temperature and Hamiltonian replica exchange and adiabatic decoupling with translational temperature scaling, were compared for a system of CaSO4 in water. It took 11 replicas in the case of temperature replica exchange to make use of a diffusion coefficient that is a factor of 1.5 larger at 350 K compared to that at 300 K. Thirty replicas were required in the Hamiltonian replica exchange with charge reduction to reach uncharged ions that have a diffusion coefficient that is 2-7 times larger than the fully charged ones. The adiabatic decoupling technique with translational temperature scaling yielded a diffusion coefficient that was 15 times larger while keeping the distribution of the water molecules around the ions unaltered with respect to the standard temperature simulation. This result illustrates the efficiency of the adiabatic decoupling technique to enhance configurational sampling.
’ INTRODUCTION The role of (bio)molecular simulation has steadily grown over the past few decades, and simulation has become a standard technique to describe and study the properties and behavior of (bio)molecular systems.1-5 When performing a simulation of a solute, polymer, or protein in aqueous solution, four aspects have to be considered:1 (i) the degrees of freedom that one chooses to be explicitly treated, (ii) the force field representing the interaction governing the motion along these degrees of freedom, (iii) the method used for an adequate sampling of those parts of configurational space that are accessible to the system at the temperature of interest, and (iv) the boundary conditions that mimic the interactions of the explicitly treated degrees of freedom with those outside the system. Here we shall focus on how to enhance sampling of configuration space for ions. The configurational sampling of a liquid is rather easy, due to the fact that it consists of many identical molecules that may exchange their position in space, and due to the relative simplicity of their interaction function. In contrast, the conformational sampling of a polymer such as a protein or a carbohydrate is much more difficult due to the connectivity of its covalent topology, which impedes a fast exchange of atom positions, and due to the complexity of its free energy surface. Sampling the configurational distribution of ions around a polymer or protein or DNA is a challenge6,7 that lies somewhere in between the two mentioned cases, as the sampling of ions is governed by longrange interactions due to their charge and slow diffusion. There is a wide range of methods available to enhance sampling.8-10 The ones usable for ions are in particular temperature replica exchange,11-13 Hamiltonian replica exchange14,15 with reducing the charge16 or changing the graining level,17,18 or adiabatic decoupling combined with increased temperature,19-21 with a r 2011 American Chemical Society
reduced interaction,22-24 or with a reduced force25 between particular sets of degrees of freedom.
’ METHODS Theory. Temperature replica exchange and Hamiltonian replica exchange with reducing the charge were investigated since both an increased temperature and a reduced charge of ions lead to an extended configurational sampling reflected by larger diffusion coefficients of the ions. The use of replica exchange with changing grain level was postponed since, to our knowledge, there is not yet a coarse grained model for water or polymers available that is thermodynamically compatible with a fine grained, all atom model for such systems, a necessary condition to obtain meaningful results from a Hamiltonian replica exchange simulation. The coarse grained simulation would mainly sample configurations that are irrelevant to the fine grained, all atom configurational ensemble. Third, adiabatic decoupling of ions and water was used together with translational temperature scaling of the ions, since this was shown to improve their diffusion, taken as a measure of extent of sampling, the most. Replica Exchange. In a typical replica exchange molecular dynamic simulation, several replicas, i.e., copies of a physical system, are simulated in parallel. These replicas differ either in temperature, for temperature replica exchange T-RE, or in their Hamiltonian, for Hamiltonian replica exchange H-RE. In Hamiltonian replica exchange, we distinguish the different Hamiltonians through a so-called coupling parameter λ, i.e., Received: November 11, 2010 Revised: January 26, 2011 Published: March 09, 2011 2931
dx.doi.org/10.1021/jp110778k | J. Phys. Chem. B 2011, 115, 2931–2936
The Journal of Physical Chemistry B
ARTICLE
Table 1. Parameters of the Nonbonded Interaction Terms in the Force Field28a ions force field parameter
Ca
S
O
OW
HW 0.41
2.0
0.54
-0.635
-0.82
31.7
99.92
47.56
51.16
q (e) C61/2 10-3 (kJ mol-1 nm6)1/2
water
C121/2 10-3 (kJ mol-1 nm12)1/2
0.7057
3.616
0.8611 (with S,O)
0
1.623
0
3.068 (with Ca) 1.841 (with OW) a
The atom partial charges are denoted by q and the Lennard-Jones parameters by C6 and C12.
H(λ) or Hn H(λn). After a chosen number of time steps, exchanges of the configurations xm and xn with adjacent λ values, λm and λn that define the Hamiltonians Hm and Hn for H-RE, or with adjacent temperatures Tm and Tn for T-RE are attempted using the detailed balance condition PS0 3 tðS0 f S00 Þ ¼ PS00 3 tðS00 f S0 Þ
ð1Þ
0
where S denotes the state before the exchange, i.e., the replica xm and the one with Tn or Hn with Tm or Hm has configuration 00 has configuration xn, and S denotes the state after the exchange, i.e., the replica with Tm or Hm has configuration xn and the one the configurawith Tn or Hn has configuration xm. PS denotes 0 tional probability of state S and0 t(S f S ) the transition probability from state S0 to state S . The relative configurational 00 probabilities of states S and S are for H-RE
and for T-RE
PS00 e-½Hm ðxn Þ þ Hn ðxm Þ=ðkB TÞ ¼ -½H ðx Þ þ H ðx Þ=ðk TÞ m m B PS0 e n n PS00 e-Hðxn Þ=ðkB Tm Þ e-Hðxm Þ=ðkB Tn Þ ¼ -Hðx Þ=ðk T Þ -Hðx Þ=ðk T Þ n B n e m B m P S0 e
ð2Þ
-Δmn
¼ minð1, e
if Δmn e 0 if Δmn g 0
Þ
T h ¼ sT T l or βh ¼ βl =sT
ð9Þ
h
The Hamiltonian for this system of N h type particles and Nl l type particles is Hðph , pl , xh , xl Þ ¼ K h ðph Þ þ V hh ðxh Þ þ V hl ðxh , xl Þ þ K l ðpl Þ þ V ll ðxl Þ
ð10Þ ð3Þ
Combining eq 1 with eq 2 or eq 3, one finds for the probability p of exchange of replicas m and n pðm T nÞ ¼ tðS0 f S00 Þ ( 1 ¼ PS00 =PS0 ¼ e-Δmn
conserving collisions becomes very small. Here we have denoted the masses of the particles constituting the Nh degrees of freedom collectively by mh and the ones of the Nl degrees of freedom collectively by ml. This decoupling is more easily achieved if Nh , Nl, which applies in the presented case taking ions as h type and water molecules as l type particles. The nondynamic equilibrium properties of a system remain unaltered by this change of mass of a subset of degrees of freedom. To enhance the sampling of the Nh h type degrees of freedom, their temperature, potential energy, or force can be scaled.25 In this work we focus on the scaling of the translational degrees of freedom of the h type degrees of freedom with a factor sT such that
ð4Þ
where Kh(ph) and Kl(pl) are the kinetic energy of the h type and the l type particles, respectively, Vhh(xh) the potential energy between h type particles, Vll(xl) the potential energy between l type particles, and Vhl(xh,xl) between h and l type particles. In the adiabatically decoupled limit, the partition function of the l type particles is Z Z l l l l ll l l hl h l Zl ðxh ; βl Þ ¼ e-β K ðp Þ dpl e-β V ðx Þ e-β V ðx , x Þ dxl ð11Þ
ð5Þ
with for H-RE
and the Hamiltonian of the h type degrees of freedom can therefore be written as H h ðph , xh ; βl Þ ¼ K h ðph Þ þ V hh ðxh Þ - ðβl Þ-1 ln Zl ðxh ; βl Þ
Δmn ¼ f½Hm ðxn Þ þ Hn ðxm Þ - ½Hm ðxm Þ þ Hn ðxn Þg=ðkB TÞ
ð12Þ
ð6Þ and for T-RE Δmn ¼ ½Hðxn Þ - Hðxm Þð1=ðkB Tn Þ - 1=ðkB Tm ÞÞ
ð7Þ
The average value of a property Q(xh) at temperature Tl, ÆQæβl, can be recalculated from ÆQæβh of Q(xh) at Th using the standard unbiasing formula ÆQ æβl ¼
Adiabatic Decoupling with Translational Temperature Scaling. A particular set of Nh degrees of freedom can be
adiabatically decoupled from the other Nl degrees of freedom of the system by increasing their mass mh ¼ sm ml
ð8Þ
with sm g 1 in such a way that the transfer of kinetic energy between the decoupled degrees of freedom through momentum
ÆQ e-ðβ
l
- βh Þ½V hh - ðβl Þ-1 ln Zl
l h l -1 Æe-ðβ - β Þ½V hh - ðβ Þ ln Zl æ
æ βh
ð13Þ
βh
where the ensemble averages are over the h type degrees of freedom and R -βl K l ðpl Þ l R -βl V ll ðxl Þ l dp e dx e l h l Z ðx ; β Þ ¼ ð14Þ l hl h þβ V ðx Þ Æe æ 2932
dx.doi.org/10.1021/jp110778k |J. Phys. Chem. B 2011, 115, 2931–2936
The Journal of Physical Chemistry B
ARTICLE
where the ensemble average in the denominator is over the l type system considered degrees of freedom.25 For the ion-water l h hh l -1 l here, the reweighting factor e-(β -β )[V -(β ) lnZ ] for the h type configurations is infinite in numerical precision, and therefore no reweighting was applied. The implementation of adiabatic decoupling in software is straightforward as it only requires scaling of atomic masses of a particular set of atoms. In addition, it should be possible to separately couple the temperature of this set of atoms to a heat bath of a given temperature. Simulated System and Simulation Procedure. CaSO4 solvated in water was chosen as a test system for comparing the sampling efficiency of the temperature and Hamiltonian replica exchange and the adiabatic decoupling with translational temperature scaling. A cubic box with an edge length of 4.082 nm was filled with 2198 SPC26 water molecules, resulting in a density of 970 kg/m3, corresponding to the density of liquid SPC water at 300 K and 1 atm.27 To this system 2 Ca2þ and 2 SO42- were added. The nonbonded interaction parameters are given in Table 1. The ideal bond length S-O in SO42- is 0.15 nm, and the O-S-O angle is 109.5 with a force constant of 520 kJ mol-1. This is according to the GROMOS force field28 53A6. Periodic boundary conditions were applied. Molecular dynamics simulations were performed with the GROMOS05 simulation software package,29 modified to incorporate adiabatic decoupling. The Table 2. Configurational and Dynamic Properties of the Ions in Aqueous Solution from Standard MD Simulations at Three Temperaturesa DCa2þ
DSO42-
T (K) (10-9 m2 s-1) (10-9 m2 s-1) ΔgCa2þOW 3 100 ΔgSO42-OW 3 100 300
1.06
1.23
0.00
0.00
325
1.17
1.31
1.82
1.39
350
1.63
1.84
4.54
3.34
T, temperature; D, diffusion coefficient, see eq 18; Δg, radial distribution difference, see eq 16.
a
Figure 1. Distribution of the potential energy per replica for the 11 replicas of the temperature replica exchange simulation, starting from 300 K (left) and going to 350 K (right) in steps of 5 K.
geometry of the water molecules was constrained by applying the SHAKE algorithm30 with a relative geometric tolerance of 10-4 on the OH bond lengths, on the intramolecular HH distance of the water molecules, and on the SO bond length of the SO42ions. The nonbonded van der Waals and electrostatic interactions were calculated using triple-range cutoff radii of 0.8/1.4 nm. The short-range interactions were calculated every time step by updating the molecular pair list for distances smaller than the first cutoff radius of 0.8 nm. For the intermediate range of distances between 0.8 and 1.4 nm, the pair list was only updated every fifth time step, and at the same time the interaction was calculated and kept unchanged between these updates. The long-range electrostatic interactions beyond the outer cutoff of 1.4 nm were represented by a reaction field31,32 with εRF = 78.5, the dielectric permittivity of water at 300 K and 1 atm. The equations of motion were integrated using the leapfrog algorithm with a time step of 2 fs. The velocities of the atoms at the beginning of the simulation were assigned from a Maxwell distribution at 300 K. For the temperature replica exchange, 11 replicas with temperatures raised in steps of 5 K from 300 to 350 K were used. The Hamiltonian replica exchange simulation comprised 30 λ values (0, 0.01, 0.025, 0.05, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.4.4, 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88, 0.915, 0.945, 0.975, 0.99, 1), where λ = 0 represents a system where the ions are fully charged and λ = 1 one with ions without charge. The charges of the ions were scaled linearly with (1 - λ). No softness was used for this transformation.33 For both RE simulations, exchanges based on a Metropolis Monte Carlo criterion were attempted every 2 ps for 1 ns. These simulations were performed at constant volume. The temperatures were weakly coupled34 to a bath with a relaxation time of 0.1 ps. The adiabatic decoupling was investigated with sm values of 1, 100, 200, 500, and 1000. The translational temperature scaling factors sT used for the ions were 1, 2, 3, and 5. These simulations were performed at constant pressure, incorporating a correction for sT in the contribution from the kinetic energy of the h type degrees of freedom to the pressure,25 by weakly coupling34 it to a bath of 1 atm with a relaxation time of 0.5 ps. The isothermal compressibility was set to 7.513 10-4 (kJ mol-1 nm-3)-1. The ions and water molecules were separately weakly coupled34 to two baths of 300 K with a relaxation time of 0.1 ps. During the runs, configurations of the system were saved every 2 ps. The various properties were taken from a 10 ns simulation that followed a 100 ps equilibration period. Additional long simulations at 325 and 350 K (with λ = 0) and λ = 0.01, 0.1,0.05, 1 (at 300 K) were performed which were used as reference simulations to which the replica exchange or adiabatically decoupled simulations could be compared. These simulations were performed at constant pressure by weakly coupling34 it to a bath of 1 atm with a relaxation time of 0.5 ps.
Figure 2. Left-hand panel: switches between the 11 replicas in the temperature replica exchange simulation as a function of time. Highlighted are the replicas that managed to go over the whole range of temperatures more than once. Right-hand panel: average acceptance ratio of switches between two adjacent replicas. 2933
dx.doi.org/10.1021/jp110778k |J. Phys. Chem. B 2011, 115, 2931–2936
The Journal of Physical Chemistry B
ARTICLE
The isothermal compressibility was set to 7.513 10-4 (kJ mol-1 nm-3)-1. They were also weakly coupled34 to a bath of the mentioned temperature with a relaxation time of 0.1 ps.
Table 3. Configurational and Dynamic Properties of the Ions in Aqueous Solution from Standard MD Simulations at Five Different Values of the Charge Reduction Factor λ Where the Charge of the Ions Is Scaled with (1 - λ)a DCa2þ (10-9 m2 s-1)
DSO42(10-9 m2 s-1)
0
1.06
1.23
0.00
0.00
0.01
1.01
1.00
16.44
12.21
0.1
1.55
1.54
41.56
26.06
0.5
2.02
1.88
54.27
31.13
1
7.85
2.64
60.54
31.67
λ
ΔgCa2þOW 3 100
ΔgSO42-OW 3 100
gRβ ðrÞ ¼
λ, charge scaling factor; D, diffusion coefficient, see eq 18; Δg, radial distribution difference, see eq 16.
a
During the runs, configurations of the system were saved every 2 ps. The various properties were taken from a 10 ns simulation that followed a 100 ps equilibration period. Analysis. The configurational distribution of the ions in water in the different systems was characterized by calculating their radial distribution functions. The diffusion coefficient was calculated as a measure of the rate of sampling. Radial Distribution Function g(r). The pair distribution function g(r) represents the probability of finding another atom at a distance r from a given atom, relative to the probability expected for a completely uniform distribution at the same density. It can be calculated by a simple histogram summation in radial shells over all molecules in the system using nðrÞ 4πr 2 ΔrF
ð15Þ
where R is the atom type Ca2þ or the S of SO42-, β the oxygen of the water, n(r) the number of atoms of type β around an atom of type R at distances between r and r þ Δr, with Δr = 0.01 nm, and F the density of the liquid.
Figure 3. Left-hand panel: switches between the 30 replicas in the charge reducing Hamiltonian replica exchange simulation as function of time. Replicas with a starting λ > 0.5 are dashed. Note that no replica manages to visit both λ = 0 and λ = 1. Right-hand panel: average acceptance ratio of switches between two adjacent replicas.
Table 4. Configurational and Dynamic Properties of the Ions in Aqueous Solution from Differently Strong (sm) Adiabatically Decoupled Simulations in Which the Temperature of the Ions Is Increased by Different Amounts (sT)a sm 1
100
200
500
1000
a
sT
DCa2þ (10-9 m2 s-1)
DSO42- (10-9 m2 s-1)
ΔgCa2þOW 3 100
ΔgSO42-OW 3 100
1
1.06
1.23
0.00
0.00
2
2.01
2.71
4.20
2.93
3
4.11
6.03
8.59
6.71
5
24.39
66.98
10.10
8.23
1
1.44
1.80
1.21
1.23
2
13.22
16.89
1.25
1.20
3
25.87
35.25
1.55
1.60
5 1
51.64 1.23
74.56 1.46
1.60 1.34
2.11 1.11
2
7.84
9.66
1.07
1.24
3
15.27
18.87
1.45
1.40
5
28.95
38.25
1.41
2.06
1
1.13
1.14
1.09
1.14
2
4.22
4.48
1.18
1.18
3
7.00
8.26
1.58
1.24
5 1
13.19 0.83
15.83 0.74
1.47 1.64
1.97 1.05
2
2.33
2.51
1.75
1.19
3
4.04
4.29
1.40
1.30
5
7.43
7.96
1.70
1.54
sm, mass scaling factor; sT, temperature scaling factor; D, diffusion coefficient, see eq 18; Δg, radial distribution difference, see eq 16. 2934
dx.doi.org/10.1021/jp110778k |J. Phys. Chem. B 2011, 115, 2931–2936
The Journal of Physical Chemistry B
ARTICLE
Figure 4. Mean-square dispacement of the Ca2þ (black, solid) and the SO42- (red, dotted) ions in aqueous solution for different mass and temperature scaling factors sm and sT in the adiabatically decoupled simulations with different ionic temperatures.
To approximate the difference between two radial distribution functions, the radial distribution difference was calculated as R rmax ref jgRβ ðrÞ - gRβ ðrÞjdr 0 R rmax ref ð16Þ ΔgRβ ¼ gRβ ðrÞ dr 0 where rmax is the cutoff for the radial distribution function, which was set to 1 nm. As reference radial distribution function gref Rβ(r), the one of the system with mass scaling factor sm=1 and temperature scaling factor sT = 1 was taken. Self-Diffusion Coefficient D. The diffusion coefficient D is obtained from the long-time limit of the mean-square displacement ð17Þ MSDðtÞ ¼ ðrðτ þ tÞ - rðτÞÞ2 τ, molecules according to the Einstein relation35 ðrðτ þ tÞ - rðτÞÞ2 τ, molecules D ¼ lim tf¥ 6t
ð18Þ
where r(t) corresponds to the position vector of the center of mass of a molecule at time t, and the averaging is performed over both time and molecules. Only the time window τ between 0 and 250 ps was used for the calculation of the diffusion constant. This avoids the inclusion of long time periods that are not independent of each other in the averaging.
’ RESULTS AND DISCUSSION The three approaches to enhance the configurational sampling of ions in aqueous solution—temperature replica exchange, Hamiltonian replica exchange with reducing the charges, and adiabatic decoupling combined with scaling of the translational temperature—were compared. Temperature Replica Exchange. Temperature replica exchange was considered since an increase in temperature from 300 to 350 K leads to an increase of the diffusion coefficient by a factor of 1.5 (Table 2). However, to span this temperature interval, at least 11 replicas are required for a sufficient number
of replica exchanges. The replicas have an exchange probability between 0.15 and 0.25. The overlap of their potential energy distributions is shown in Figure 1. For 1 ns of replica exchange simulation with exchange attempts every 2 ps, only 3 replicas manage to populate 300 and 350 K (bold lines in Figure 2). The gain in diffusion by including higher temperature replicas is lost by the necessity to have a sufficient number of replicas. This problem is even more pronounced for bigger systems since the number of replicas necessary is proportional to the square root of the system size.36 Simply simulating the system at a higher temperature is not an option, since the structure (see the radial distribution difference in Table 2) changes significantly if the temperature of the whole system is increased. Hamiltonian Replica Exchange Using Charge Reduction. As a reduction of the charge of the ions leads to a considerable increase in the diffusion coefficient by a factor of 2-7 (Table 3), a Hamiltonian replica exchange which reduces the ionic charge was performed. To obtain an exchange probability between 0.07 and 0.4, 30 replicas were required. In the simulated time of 1 ns, none of the 30 replicas manages to cross the charge range completely (Figure 3). As in the case of the temperature replica exchange, the gain in diffusion from replicas with λ values close to 1 is lost by the necessity to have a sufficient number of replicas. Since the radial distribution differences (Table 3) are also significant for a simulation with reduced ionic charge, a simulation with a reduced ionic charge does not offer a viable alternative to replica exchange. Adiabatic Decoupling with Translational Temperature Scaling. To get an impression how much the structure of the water molecules around the adiabatically decoupled ions that have a higher temperature has changed, the radial distribution difference ΔgRβ is shown in Table 4. Since scaling the mass with sm without scaling the temperature should not influence the radial distribution function, the results for sT = 1 can be taken as a reference value for contributions of the noise in the radial distribution function. The maximum ΔgCaOW is 1.64 observed for sm = 1000, and the maximum for ΔgSOW is 1.23 found for sm = 100. These deviations from zero are due to limited sampling and 2935
dx.doi.org/10.1021/jp110778k |J. Phys. Chem. B 2011, 115, 2931–2936
The Journal of Physical Chemistry B the finite resolution of 0.01 nm used in the representation of ΔgRβ(r). For comparison, the corresponding values when raising the temperature to 325 K are larger, 1.82 for ΔgCaOW and 1.39 for ΔgSOW. For sm = 1, i.e., without adiabatic decoupling, ΔgRβ grows fast with increasing sT. All systems with sm > 1, except two, sm = 1000 with sT = 2 or 5, display ΔgCaOW values within the variation of sT = 1. ΔgSOW seems to be a bit more sensitive to raising the temperature: all systems with sm > 1 using sT = 3 or 5 yield ΔgSOW values that are somewhat larger than the variation 1.23 found for sT = 1. However, ΔgSOW for sT = 5 nicely decreases for increasing sm, i.e., for better adiabatic decoupling. The diffusion coefficients in Table 4 show a remarkable increase for the adiabatically decoupled systems, sm > 1, with an increased temperature, sT > 1. The diffusion increases with increasing sT and decreases with increasing sm. This implies that an optimum is to be found between sufficient decoupling and not too large ionic mass. For sm > 1 and sT > 1, the time dependence of the mean-square displacement (Figure 4) is not linear anymore. This implies that the values for D in Table 4 would change if the considered time window for obtaining the D values from the mean-square displacement, eq 18, would be varied. Thus the D values are mainly given for ease of comparison.
’ CONCLUSIONS Adiabatic decoupling with translational temperature scaling offers a method to considerably increase the diffusion of particular degrees of freedom, e.g., ions in aqueous solution, without perturbing the configurational distribution of the surrounding degrees of freedom, i.e., water, significantly. An increase of a factor of 15 in the diffusion coefficient was observed with a mass scaling factor sm = 100 and a temperature scaling factor sT = 2, which represents a system for which almost no difference in the radial distribution function compared to the system at sm = 1, sT = 1 was observed. To achieve this, only one simulation had to be performed compared to the 11 replicas required in temperature replica exchange and 30 replicas in Hamiltonian replica exchange with charge reduction. Adiabatic decoupling with temperature scaling of the ionic degrees of freedom is, therefore, a very promising method to enhance the configurational sampling of ions around a solute, such as a charged protein, carbohydrate, RNA, or DNA, in aqueous solution. ’ AUTHOR INFORMATION Corresponding Author
*E-mail
[email protected].
’ ACKNOWLEDGMENT The authors thank Denise Steiner for help with the replica exchange simulations. This work was financially supported by the National Center of Competence in Research (NCCR) in Structural Biology, by grant number 200020-121913 of the Swiss National Science Foundation, and by grant number 228076 of the European Research Council, and by a gift of Sika group to the ETH foundation, which is gratefully acknowledged. ’ REFERENCES
ARTICLE
(2) Scheraga, H. A.; Khalili, M.; Liwo, A. Annu. Rev. Phys. Chem. 2007, 58, 57–83. (3) van Gunsteren, W. F.; Dolenc, J. Biochem. Soc. Trans. 2008, 36, 11–15. (4) Lindahl, E.; Sansom, M. S. P. Curr. Opin. Struct. Biol. 2008, 18, 425–431. (5) Khalili-Araghi, F.; Gumbart, J.; Wen, P.-C.; Sotomayor, M.; Tajkhorshid, E.; Schulten, K. Curr. Opin. Struct. Biol. 2009, 19, 128–137. (6) Peric, L.; Pereira, C. S.; Perez, S.; H€unenberger, P. H. Mol. Simul. 2008, 34, 421–446. (7) Peric-Hassler, L.; H€unenberger, P. H. Mol. Simul. 2010, 36, 778–795. (8) Christen, M.; van Gunsteren, W. F. J. Comput. Chem. 2008, 29, 157–166. (9) Liwo, A.; Czaplewski, C.; Oldziej, S.; Scheraga, H. A. Curr. Opin. Struct. Biol. 2008, 18, 134–139. (10) Christ, C. D.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2010, 31, 1569–1582. (11) Swendsen, R. H.; Wang, J. S. Phys. Rev. Lett. 1986, 57, 2607–2609. (12) Hukushima, K.; Nemoto, K. J. Phys. Soc. Jpn. 1996, 65, 1604–1608. (13) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141–151. (14) Sugita, Y.; Kitao, A.; Okamoto, Y. J. Chem. Phys. 2000, 113, 6042–6051. (15) Woods, C. J.; Essex, J. W.; King, M. A. J. Phys. Chem. B 2003, 107, 13703–13710. (16) Christ, C. D.; van Gunsteren, W. F. J. Chem. Phys. 2007, 126, 184110. (17) Lwin, T. Z.; Luo, R. J. Chem. Phys. 2005, 123, 194904. (18) Christen, M.; van Gunsteren, W. F. J. Chem. Phys. 2006, 124, 154106. (19) Rosso, L.; Tuckerman, M. E. Mol. Simul. 2002, 28, 91–112. (20) VandeVondele, J.; R€othlisberger, U. J. Phys. Chem. B 2002, 106, 203–208. (21) Rosso, L.; Minary, P.; Zhu, Z. W.; Tuckerman, M. E. J. Chem. Phys. 2002, 116, 4389–4402. (22) Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13749–13754. (23) Liu, P.; Huang, X. H.; Zhou, R. H.; Berne, B. J. J. Phys. Chem. B 2006, 110, 19018–19022. (24) Huang, X. H.; Hagen, M.; Kim, B.; Friesner, R. A.; Zhou, R. H.; Berne, B. J. J. Phys. Chem. B 2007, 111, 5405–5410. (25) Kunz, A.-P. E.; Liu, H.; van Gunsteren, W. F. J. Chem. Phys. 2010submitted for publication. (26) 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.; Reidel: Dordrecht, Germany, 1981. (27) Gl€attli, A.; Daura, X.; van Gunsteren, W. F. J. Chem. Phys. 2002, 116, 9811–9828. (28) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2004, 25, 1656–1676. (29) Christen, M.; H€unenberger, P. H.; Bakowies, D.; Baron, R.; B€urgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenholz, M. A.; Kr€autler, V.; Oostenbrink, C.; Peter, C.; Trzesniak, D.; van Gunsteren, W. F. J. Comput. Chem. 2005, 26, 1719–1751. (30) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (31) Barker, J. A.; Watts, R. O. Mol. Phys. 1973, 26, 789–792. (32) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. J. Chem. Phys. 1995, 102, 5451–5459. (33) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222, 529–539. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (35) Berendsen, H. J. C. Simulating the physical world, 1st ed.; Cambridge University Press: Cambridge, U.K., 2007.
(1) van Gunsteren, W. F. Angew. Chem., Int. Ed. 2006, 45, 4064–4092. 2936
dx.doi.org/10.1021/jp110778k |J. Phys. Chem. B 2011, 115, 2931–2936
