A pH-Replica Exchange Molecular Dynamics method is proposed to improve the coupling between the conformational and protonation samplings. By swapping the conformations between two neighbor replicas, which are at different pHs, this method increase the conformational sampling directly and protonation state sampling indirectly. The method has been applied on seven biological systems. The pHEX MD method, predicted the pKa values of small model compound correctly, and converged faster than Constant pH MD method. Further, pHEX has been applied on ADFDA polypeptide (without capping) and a hepta-peptide derived from the ovomucoid third domain (OMTKY3). In all tested compounds, the predicted pKa by the pHEX is very close to CpH MD, however the pHEX has the advantage of more efficient conformational and protonation state samplings.
Introduction
Function, dynamics and structure of the proteins are highly correlated with Solution pH. Many important biological phenomena function only at a certain range of pH, the examples are include protein folding/misfolding1-3, enzyme catalysis4-6 and ligand-substrate docking7,8. The pH change shifts the protonation equilibrium of titrable residues, which leads to the change in the conformation of the protein itself.9-12
The pH value, at which, the ratio of deprotonated state of a titrable residue (usually side chain) to protonated state of that becomes one, is the pKa of that titrable residue. The pKa of an ionizable side chain in a protein is highly dependent on the electrostatic environment of the side chain, which is affected by the prtotonaion states of other titrable groups and the conformation state of the protein. Because the change in protonation state is involving bond breaking, which is a quantum phenomenon, essentially in traditional Molecular Dynamics (MD) methods, protonation state of ionizable species are constant. The constant protonation MD methods suffer from two big disadvantages.13 First, at pH solution near the pKa of any of the titrable residues, the protonation state of that residue is not constant and this decouples between sampling of protonation and conformation states. Second, fixing the protonation state of a titrable residue required to know the pKa of its ionizabe side chain, which is priori unknown.
Due to distinction of pH solution as an external and controllable variable in experimental methods, like pressure, temperature, etc., the importance of constant pH MD has been recognized. In the last two decades, many Constant pH MD methods have been developed.14,15 The goal of CpH MD is to describe correctly the protontion equilibrium coupled by conformation equilibrium at certain pH. The majority of CpH MD methods can be divided to continuous and discrete protonation state methods.
Continuous protonation methods use continuous protonation parameter to peturbe ionizable residue from protonated to deprotonated state. In 1994, Mertz and Pettitt16 developed a grand canonical method for simulating a simple chemical reaction. They applied the method for exchanging the proton between water molecules and an ionizable side chain. In 1997, Baptista et al.14 introduced continuous constant pH method in implicit solvent base on meanfield approximation. In 2001, Börjesson et al.15 used a weakly coupled proton bath to continuously adjust the protnation fraction of each titrable group toward the equilibrium. More recently, Brooks group has developed the continuous protonation state method further.17-22 In the case of highly coupled titration group, where cooperativity effect is non-negligible, this model leads to inappropriate estimation of physical variables. To alleviate this problem, Lee et al.22 has introduced λ-dynamics23 along the protonation coordinate, where λ equal to zero and λ equal to 1 correspond to fully protonated and deprotonated states of a titrable group respectively. They have added a biasing potential, which is a function of λ and centered at λ equal 0.5, to favor protonation state coordinate value to fully protonated/deprotonated states.
In contrast to continuous protonation state method, in discrete models, the protonation state of the ionizable group is either zero or one during the simulation.
These models use a hybrid MD-MC scheme; while the MD is used to sample the conformation space, At every certain MD steps, A Metropolis MC24 question will be asked for changing the protonation state/states. Many variation of discrete protonation state have been developed.13,25-34 While Baptista group 26-29,35used the explicit solvent for propagtion of coordinates and Poisson-Boltzmann (PB) method for calculation of energy in MC parts, Walczak et al.34 employed the Lengevin Dynamics for MD parts and the PB method for MC steps. The big disadvantage of these approaches are that the energy calculation methods are different in protonation and conformation sampling parts. For tackling this problem, Bürgi et al.30 applied the Thermodynamics Integration method (TI) to calculate the transition energy between the prtotonated and deprotonated species in explicit solvent MD. There are at least two drawbacks for this approach. First, the computational cost of TI calculation is demanding, and this limit the number of MC sampling. Second, even in the case of MC step rejection TI calculation will perturb the MD part. In 2004, Mongan et al. developed a discrete protonation MD method by applying the GB implicit solvent method in both MD and MC sampling part.13 A more descriptive explanation of this method can be found in method section. This method is implemented in AMBER Molecular Dynamics Package.36
Rather than protonation space, efficient and accurate sampling of conformation space of proteins always been challenging. Many theoretical methods proposed to overcome the free energy barriers in conformational space. Among those ,generalized ensemble methods37-40 such as Replica Exchange MD (REMD)41, Simulated Tempering42 (ST) and multicanonical ensemble43,44, are well stablished. In those methods, in addition to conformational space (i.e., regular MD) the system perform a random walk in energy or temperture space for reducing the kinetic trapping problem. REMD has at least two advantages respect to the other methods. First, its algorithm is highly parallelizable. Second, its weighting factor, e.g. Boltzmann factor, is priori known.
Due to the coupling of protonation and conformation sampling, recently, it has been tried to combine the Constant pH MD with REMD method. Khandogin et al.18 combined the continuous protonation Constant pH method with REMD algorithm. They applied Their method, which has been named REX-CPHMD, for pKa prediction protein folding, pH dependent conformation, etc. Recently, Meng and Roitberg45 utilized a hybrid method by combining the Temperature REMD (T-REMD) and discrete protonation Constant pH MD.
Here for the first time, we intoduce pH-Exchange MD by combining the discrete protonation Constant pH MD(proposed by Mongan et al.13) and Hamiltonian REMD (H-REMD)[ref].
We tested our method by applying it on five depeptides, sequence of Ala-Asp-Phe-Asp-Ala (ADFDA), and a heptapeptide derived from OMTKY3.Because two ends of ADFDA are not capped, two Asp residues have different electrostatic environments and essentially their pKas deviate from the pKa of unperturbed of Asp sidechain. Previously Dlugosz and Antosiewicz31,32 calculated the pKa of the heptapeptide derived from ovomucoid third domain (OMKTY3). Our goal is to show that the pHEX MD improve the sampling efficiency in both protonation and conformation spaces.
Hereafter, Constant pH MD (CpH) is referring to the Mongan's CpH MD approach, unless, it has been mentioned .
Method
Constant pH molecular dynamics
The goal of CpH MD is sampling the equilibrium between protonated and deprotonated state of titrable sites at a given pH. The free energy difference between protonated and deprotonated states rules the ratio of concentration of them. This free energy difference can not be calculated by Molecular Mechanics (MM), because, Changing from protonated state to deprotonated state is involving bond breaking/forming phenomena, which required quantum calculation. Mongan et al.13,46 proposed a method, which uses the pre-calculated pKa of reference compounds, to overcome this barrier. Reference/model compounds are AMBER terminology, which uses for depeptide contain titrable residues (i.e., Ace-titrable residue-Nme). In the paper, Mongan assumed the free energy difference between two forms of protonation states can be divided into QM and MM parts.
Equation 1
Further, they assumed QM component is the same for titrable site in protein and reference compound i.e.,
.
Equation
Since pKa of reference compound is known, the difference in free energy of its deprotonated and protonated forms at a given pH is
.
Equation
From Equation and Equation QM part of difference in free energy of change in protonation state of reference model can be written as,
.
Equation
Using Equation and Equation , the transition free energy of titrable residue from deprotonation state to protonation state can be written as,
,
Equation
Where ï„Gprotein,MM is the molecular mechanics part of free energy of the titrable site of the protein. Using equation 5, there is no need to calculating the QM contribution of free energy in protein. This method is adopted in AMBER MD suite36, but now, it is only available in GB implicit solvent model. ï„Gprotein,MM is the difference between the potential energies of the protein titrable residue at the current protonation state and new protonation state. pKa,ref and ï„Gref,MM are pre-calculated parameters in AMBER. Every a few MD steps a Metropolis MC question24 will be asked to change the protonation state of the titrable residue and ï„Gprotein will be used to decide about accepting or rejecting MC move. In other word, MC move take care of sampling in protonation space, and MD steps take care of sampling in configuration space. During MD steps protonation state is constant.
Titration curve:
The pKa of a titrable residue is related to pH environment through the Henderson-Hasselbalch (HH) equation,
,
Equation
Here and are deprotonated and protonated concentrations respectively; n is the Hill coefficient. In the case of non-interacting ionizabe residue n is equal to 1; In the case of interacting titrable residue it deviate from one because of cooperativity47.
Because of ergodicity assumption in MD, the ratio of the time that titrable residue spend at deprotonated states to the time that of it spend in protonated state, can be considered as a ratio of concentration of deprotonated state to that of protonation state.
pH- Exchange(pHEX):
In contrast to regular temperature replica exchange41, which each replica runs in different temperature, in pHEX, each replica runs in different pH but the same temperature.
Each replica run constant pH MD at a unique pH, and periodically an exchange of conformation between two adjacent replicas is attempted.
According to the equilibrium condition of states, the detailed balance condition can be written as
.
Equation
At the exchange moment, if we only swap the conformations between two replica and keep the protonation state unchanged, then the generalized states of X and X' can be written as
and .
Here and are the conformation and protonation states respectively and the ith column is related to the ith replica .
Because all N replicas are independent, the probability of the system being in generalized state of can be written as , where is the probability of replica i at conformation and protonation state of .
Substituting those in Equation will result
.
Equation
In the canonical ensemble regime (NVT) this can be written as
,
Equation
where and ï¢=KBT.
Then the exchange probability according to Metropolis MC criteria is
Because the two exchanging replicas are runing at the same temperature, in computing of ∆, the kinetic energy terms will cancel each other, so only the potential energies are required.
Another possibility is trying to change both the conformation and the protonation states between two replicas. But in this case MC question will be trivial, because ï„ is equal to zero, which means every exchange will be accepted.
Simulation details
For testing our method, we choose and ran simulations on three categories of systems. First, titrable model compounds (Ace-titrable residue-Nme), those used as reference compounds (see the method part) in the CpH method in AMBER 1036.
Simulation times, for all model compounds, were three ns in both CpH MD and pHEX MD (for each replica) methods. We used eight replicas for all the model compounds in pHEX.
We also tested our simulation method on a pentapeptide model, Ala-Asp-Phe-Asp- Ala (ADFDA), which has no cap at its two ends. Because two ends of ADFDA are not capped, two Asp experience different electrostatic environments and have different pKas. The simulation times of ADFDA were ninety ns and ten ns for CpH and pHEX (for each replica) methods respectively. Eight replicas, at pH=2.5-6.0 with increment of 0.5, have been used in pHEX method.
The third system heptapeptide derived from OMTKY3 (ACE-Ser-Asp-Asn-Lys-Thr-Tyr-Gly-NME). Dlugosz and Antosiewicz31,32 studied this heptapeptide and they predicted the pKa of 4.24 for Asp3 using their CpH MD method. The simulation time of 100 ns and 10 ns used for CpH and pHEX methods respectively. Twelve replicas, at pH=2-13 with increment of 1.0, have been used in pHEX method.
All calculations in this paper are done using AMBER 10 molecular simulation suite. The AMBER ff99SB force field48 and OBC Generalized Born implicit solvent49 (igb=2) have been used in all simulations. Shake algorithm50 applied to constrain the connected bonds to hydrogen atoms, which allowed the use two ps MD steps. The cutoff for non-bonded interaction and the born radii were chosen 30 and 999 Å respectively in all calculations.
In all pHEX simulations, except the heptapeptide derived from OMTKY3, replicas have been allowed to ask for the conformation exchange every 1000 MD steps. In the case of heptapeptide, for having a faster structural convergence, the conformational exchanges attempted in every 500 MD step. The exchange ratios between replicas for all system were between 0.35 and 1.
Cluster analysis and Local convergence
In the case of ADFDA, Cluster analysis was done using Dr. Simmerling group's program [ref] to test the structural convergence. For all the compounds the local conformational sampling was studied by comparing the probability density of the backbone dihedral angles (ïª and ï¹) for all the residues. Probability densities were computed by dividing (ïª,ï¹) 2d-space to 10ï‚°ï‚´ 10ï‚° bins. The resulted histograms were normalized and contours were plotted.
Discussion
Titrable model compounds
We applied our pH-exchange method to titrable model compounds. The simulation time has been chosen 3 ns for both pHEX (for each replica) and CpH methods. As Table shows there are good agreements among experimental, CpH and pHEX methods for model compounds pKa. All pKa values have been calculated by Henderson-Hasselbalch equation and linear regression technique.
Table .pKa of Model compounds calculated by different methods
ASP
GLU
HIS
TYR
LYS
Experimental value
4.0
4.4
6.3
9.6
10.5
CpH MD
3.9
4.4
6.3
9.6
10.4
pHEX MD
3.9
4.4
6.4
9.7
10.4
As an example the Hill plots of caped Lysine has been plotted for both pHEX and CpH methods, which show a great conformity.
Figure . Comparison of titration curve between pHEX and CpH methods for Lys reference model. The Asymptotic Standard Errors for slopes of the fitted lines are ±0.014 and ±0.020 for CpH and pHEX methods respectively.
ADFDA model compound
The Second test applied on a model peptide, Ala-Asp-Phe-Asp- Ala (ADFDA), which its two ends are not capped. Because of that, two Asp ionizable side chains have different electrostatics environments. Asp2 is close to the NH3+, which means at pH equal 4.0, its deprotonated state is more favored and its pKa is shifted below 4.0; Asp4 is close to COO-, which means its deprotonated state at pH equal 4.0 is less favored and its pKa is shifted above 4.0.
Figure . The titration curves of Asp side chains in ADFDA calculated by both constant pH MD (blue and purple) and pH-Exchange (red and green); The titration curve of Asp reference compound (light blue).
The titration curves of Asp2 and Asp4 have been plotted. As Figure shows the titration curves of both Asp2 and Asp4 show a shift to the left and right respectively compare to the Asp model compound. In the case of Asp4, there is a good agreement between pHEX and CpH methods, however in the case of Asp 2 the difference between the titration curves of two methods is more. The pKa and Hill coefficent of Asp2 and Asp4 for the both methods have been calculated by linear fitting on Hill plot. As it proposed by Table , the results of the both methods are close.
Table . pKa prediction and Hill coefficient of fitted from the HH equation.
ASP2
ASP4
pKa
Hill coefficient
pKa
pHEx
3.77
1.02
4.47
CpH
3.75
0.89
4.47
As an example, in Figure , the cumulative protonation fraction of Asp2 and Asp4 for both methods at pH=4.0 has been plotted. According to Figure , for the both titrable groups pHEX shows a faster convergence in protonation space respect to regular CpH method.
Figure . Cumulative protonation fractions of Asp side chain in ADFDA VS MC titration steps at pH=4.0.
Figure . Left: Cluster population of CpH VS pHEX runs at pH=4.0. Right: Cluster population of 2 different runs pHEX at pH=4.0.
Global conformation convergence of pHEX has been studied by cluster analysis method. Figure , on right, shows the comparison between the cluster populations of two 10 ns runs of pHEX simulation at pH=4.0, which have used two different random seeds. As correlation coefficient between two pHEX runs (R2 value) proposes both simulation has converged. In Figure , on left, the percentages of populations of CpH at pH=4.0 has been plotted to that of pHEX method. The correlation coefficient between CpH and pHEX methods for the rest of pHs are listed in Table .Although the length of simulation times in CpH method and pHEX are 90 and 10 ns respectively, but the correlation coefficient values indicate that both methods sampled the same conformational space and generated the same structure ensemble.
Table . Correlation Coefficient between percentage of cluster populations of CpH and pHEX methods.
pH
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
Correlation Coefficient (R2 )
0.96
0.92
0.94
0.96
0.96
0.93
0.98
0.96
Figure . Ramachandran plots for Asp4 in ADFDA. Left: CpH method. Right: pHEX method
For having better resolution in convergence, local convergence of dihedral angle has been considered. As Figure shows, both methods sampled the same conformation in Φ and Ψ angles of backbone of Asp4.
Heptapeptide derived from OMTKY3
We also applied pH-Exchange method to a heptapeptide derived from OMTKY3 (ACE-Ser-Asp-Asn-Lys-Thr-Tyr-Gly-NME).
According to NMR experiments33,51, the pKa of Asp in this heptapeptide is about 3.6, which shows a -0.4 shift in pKa unit respect to that of Asp model compound.
Dlugosz and Antosiewicz's constant pH MD method predicted the pKa shift of + 0.24 for this heptapeptide, while we measured the pKa of 3.7 and 3.6 using the discrete pH MD and pH-Exchange methods respectively. Those pKa predictions are in excellent agreement with experimental value.
Rather than Asp3, there are two more titrable groups in the heptapeptide, Lys5 and Tyr7. As Figure and Figure propose, the titration curves of Asp3, Lys5 and Tyr7 are closely follow the same trend in both pHEX and CpH methods.
Figure . The titration curves of Asp3 in the heptapeptide derived from OMTKY3.
Figure .The titration curves of Lys5 and Tyr7 in the heptapeptide derived from OMTKY3.
Figure . Hill plot of Tyr7.
Applying HH equation and the linear regression method, the pKa of the three titrable groups has been calculated. Figure shows the comparison between Hill plots of CpH and pHEX methods for Tyr 7. The great conformity of two plots is clear.
The pKa predictions of Constant pH MD are 3.7, 10.6 and 9.9 for Asp3, Lys5 and Tyr7 respectively, while those of pH-Exchange method are 3.6, 10.1 and 10.6.
Table . Comparison between pHEX and CpH for pKa values of titrable residues in the heptapeptide derived from OMKTY3
Asp3
Lys5
Tyr7
pHEx
3.6
10.8
10.1
CpH
3.7
10.6
9.9
Figure .cumulative protonation state VS MC titration steps (pH=8.0,9.0,10.0,11.0,12.0 and 13.0) for TYR7 .The comparison between CpH(left ) and pHEX(right ).
To compare the convergence speed in protonation/deprotonation spaces of titrable groups in both methods, we studied the cumulative protonation/deprotonation fraction as function of MC titration steps for all ionizable residues. As an example, in the case of Tyr 7, Figure shows in protonation space, pHEX MD converges faster than CpH MD.
Figure .The difference between instant Cumulative Root Mean Square Deviation of structures (C-RMSD(t)) and C-RMSD(10 ns) VS simulation time for both CpH and pHEX methods at pH=4.0 (Left) and pH=10 (Right).
The conformational convergence has been studied by calculation of instant Cumulative Root Mean Square Deviation of structures from the first snapshot (C-RMSD(t)) for CpH and pHEX methods. Because of frequent change in the protonation states of titrable groups at pH=4.0 and 10.0, here, the C-RMSD(t) of have been presented for both method at those pHs to show the performance of two methods. In both plots, C-RMSD (t=10 ns) has been subtracted from C-RMSD(t) for better visualization purposes. In both pHs, pHEX MD shows a faster and smoother conformational convergence.
Figure . Ramachandran plots for Asp3 in the hepta-peptide at pH=4.0. Left: CpH method. Right: pHEX method
The quality of local structure sampling of two methods studied by considering the Ramachandran plots for all residues. In Figure , as an example, the Ramachandran plots of CpH (Left) and pHEX (Right) for Asp3 at pH=4.0 have been presented. Comparing two plots shows clearly that pHEX methods sampled more area in two-dimensional space of ï¹-ï¦. After 100 ns of simulation in CpH, it hardly sampled the Beta-Sheet area, while after only 10 ns of pHEX simulation, it not only sampled better Beta-sheet area, but also it touched left handed alpha helix part, which never touched by CpH MD method.
The sampling results propose that pHEX accelrate both conformational and protonationan convergence for the heptapeptide.
Conclusion
In presented work, Hamiltonian Replica Exchange Molecular Dynamics has been combined with Constant pH MD with discrete protonation state model. The method applied on seven small peptides to compare with CpH MD. The predicted pKa by pHEX MD is in great agreement with experimental and CpH methods. Compare to CpH MD, pHEX MD converges faster in both conformational and protonation spaces. In contrast to the Temperature Replica Exchange Molecular Dynamics methods (T-REMD), in pHEX the replica ladder (pH) is very limited, so the overlap of energy and consequently the exchange ratio between neighbor replicas are always high. It is expected in the systems with highly coupled confomational and protonation states, like proteins with pH-dependent structures, sampling efficiency of this method increases. Base on this, our future is included the application of pHEX MD on the systems with highly pH-dependent dynamics.
1. Bierzynski, A., Kim, P.S. & Baldwin, R.L.Proc. Natl. Acad. Sci. U. S. A. 79, 2470 (1982).
2. Shoemaker, K.R. et al.Proc. Natl. Acad. Sci. U. S. A. 82, 2349 (1985).
3. Schaefer, M., Van Vlijmen, H.W.T. & Karplus, M.Adv. Protein Chem. 51, 1 (1998).
4. Demchuk, E., Genick, U.K., Woo, T.T., Getzoff, E.D. & Bashford, D.Biochemistry 39, 1100 (2000).
5. Dillet, V., Dyson, H.J. & Bashford, D.Biochemistry 37, 10298 (1998).
6. Harris, T.K. & Turner, G.J.IUBMB Life 53, 85 (2002).
7. Antosiewicz, J., Briggs, J.M. & McCammon, J.A.Eur. Biophys. J. Biophy. 24, 137 (1996).
8. Hunenberger, P.H., Helms, V., Narayana, N., Taylor, S.S. & McCammon, J.A.Biochemistry 38, 2358 (1999).
9. Hill, T.L.J. Am. Chem. Soc. 78, 3330 (1956).
10. Simonson, T., Carlsson, J. & Case, D.A.J. Am. Chem. Soc. 126, 4167 (2004).
11. Tanford, C. & Kirkwood, J.G.J. Am. Chem. Soc. 79, 5333 (1957).
12. Warshel, A.Nature 330, 15 (1987).
13. Mongan, J., Case, D.A. & McCammon, J.A.J. Comput. Chem. 25, 2038 (2004).
14. Baptista, A.M., Martel, P.J. & Petersen, S.B.Proteins 27, 523 (1997).
15. Borjesson, U. & Hunenberger, P.H.J. Chem. Phys. 114, 9706 (2001).
16. Mertz, J.E. & Pettitt, B.M.Int. J. Supercomput. Ap. 8, 47 (1994).
17. Khandogin, J. & Brooks, C.L.Biophys. J. 89, 141 (2005).
18. Khandogin, J. & Brooks, C.L.Proc. Natl. Acad. Sci. U. S. A. 104, 16880 (2007).
19. Khandogin, J. & Brooks, C.L.Biochemistry 45, 9363 (2006).
20. Khandogin, J., Chen, J.H. & Brooks, C.L.Proc. Natl. Acad. Sci. U. S. A. 103, 18546 (2006).
21. Khandogin, J., Raleigh, D.P. & Brooks, C.L.J. Am. Chem. Soc. 129, 3056 (2007).
22. Lee, M.S., Salsbury, F.R. & Brooks, C.L.Proteins: Struct., Funct., Bioinf. 56, 738 (2004).
23. Kong, X.J. & Brooks, C.L.J. Chem. Phys. 105, 2414 (1996).
24. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. & Teller, E.J. Chem. Phys. 21, 1087 (1953).
25. Baptista, A.M.J. Chem. Phys. 116, 7766 (2002).
26. Machuqueiro, M. & Baptista, A.M.J. Phys. Chem. B 110, 2927 (2006).
27. Machuqueiro, M. & Baptista, A.M.Biophys. J. 92, 1836 (2007).
28. Machuqueiro, M. & Baptista, A.M.Proteins: Struct., Funct., Bioinf. 72, 289 (2008).
29. Machuqueiro, M. & Baptista, A.M.J. Am. Chem. Soc. 131, 12586 (2009).
30. Burgi, R., Kollman, P.A. & van Gunsteren, W.F.Proteins 47, 469 (2002).
31. Dlugosz, M. & Antosiewicz, J.M.Chem. Phys. 302, 161 (2004).
32. Dlugosz, M. & Antosiewicz, J.M.J. Phys. Chem. B 109, 13777 (2005).
33. Dlugosz, M., Antosiewicz, J.M. & Robertson, A.D.Phys. Rev. E 69, 021915 (2004).
34. Walczak, A.M. & Antosiewicz, J.M.Phys. Rev. E 66, 051911 (2002).
35. Baptista, A.M., Teixeira, V.H. & Soares, C.M.J. Chem. Phys. 117, 4184 (2002).
36. Case, D.A. et al. AMBER 10. null, (2008).
37. Li, H.Z., Fajer, M. & Yang, W.J. Chem. Phys. 126, 024106 (2007).
38. Mitsutake, A., Sugita, Y. & Okamoto, Y.Biopolymers 60, 96 (2001).
39. Zheng, L.Q., Chen, M.G. & Yang, W.Proc. Natl. Acad. Sci. U. S. A. 105, 20227 (2008).
40. Zheng, L.Q., Chen, M.G. & Yang, W.J. Chem. Phys. 130, 234105 (2009).
41. Sugita, Y. & Okamoto, Y.Chem. Phys. Lett. 314, 141 (1999).
42. Lyubartsev, A.P., Martsinovski, A.A., Shevkunov, S.V. & Vorontsovvelyaminov, P.N.J. Chem. Phys. 96, 1776 (1992).
43. Berg, B.A. & Neuhaus, T.Phys. Lett. B 267, 249 (1991).
44. Berg, B.A. & Neuhaus, T.Phys. Rev. Lett. 68, 9 (1992).
45. Meng, Y. & Roitberg, A.E. Constant pH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model. Journal of Chemical Theory and Computation 6, 1401-1412 (2010).
46. Mongan, J. & Case, D.A.Curr. Opin. Struct. Biol. 15, 157 (2005).
47. Stephens, E., O'Brien, D.P., Zhou, M. & Williams, D.H. Understanding noncovalent interactions: ligand binding energy and catalytic efficiency from ligand-induced reductions in motion within receptors and enzymes. Angewandte Chemie International Edition 43, 6596-616
48. Hornak, V. et al.Proteins: Struct., Funct., Bioinf. 65, 712 (2006).
49. Ryckaert, J.P., Ciccotti, G. & Berendsen, H.J.C.J. Comput. Phys. 23, 327 (1977).
50. Onufriev, A., Bashford, D. & Case, D.A.J. Phys. Chem. B 104, 3712 (2000).
51. Dlugosz, M. & Antosiewicz, J.M.J. Phys.: Condens. Matter 17, S1607 (2005).
