Constant pH molecular dynamics study on the doubly mutated staphylococcal nuclease: capturing the microenvironment

Abstract

A constant pH molecular dynamics (CpHMD) simulation method utilizing generalized Born electrostatics is employed to examine conformational changes of a protein containing a buried ionizable pair in its hydrophobic core, specifically considering a hyperstable variant of staphylococcal nuclease (Δ+PHS). It is observed that the average geometrical parameters obtained from CpHMD simulation at pH 8 and from quantum mechanical calculation nicely match to the reported crystal structure of V23E/L36K variant of Δ+PHS. Strong hydrogen bond formation and salt bridge interactions in the microenvironment around the Glu-23/Lys-36 pair stabilize the core of the V23E/L36K variant. Small rearrangements of the backbone carbonyl oxygen of Gly-20 and the side chain of hydroxyl group of Thr-62 are observed, and these changes in orientation allow for interaction with Lys-36 and Glu-23, respectively, in V23E/L36K variant. We do observe a long hydrogen bond correlation and lifetime (716 ps) in the microenvironment. All these observations suggest that the V23E/L36K variant is stable and CpHMD method can capture small rearrangements of the residues in a protein.

---

Introduction

Recently, with the advent of sophisticated experimental techniques and advanced molecular dynamics (MD) protocol, the folding process of a protein can be probed at the molecular level by monitoring many molecular level interactions, namely, van der Waals interactions, hydrogen bonds, disulfide bridges, salt bridges, hydrophobic interactions, etc.^1–15 It is well known that the conformational rearrangements of a protein are often associated with their functions. The conformational rearrangements of a protein are influenced by the external parameters such as temperature, solvent, solvent pH, salt concentration, presence of other molecules, etc. The folded state of a particular protein is preferred by minimizing the contact of neutral hydrophobic residues towards the water and maximizing the exposure of charged residues towards the solvent,^1,2 thus making the protein core quite dry. Significant difference is also observed in polarizability between the confined water inside a folded protein and the bulk water, and giving the protein dielectric heterogeneity.^15–18 Recent theoretical study also shows significant difference of dielectric constant between confined water cluster and bulk water.¹⁹

A nonpolar residue can easily be mutated with a polar charged residue in the dry core of a protein, namely, staphylococcal nuclease (SNase) without suffering any significant changes to the global structure and stability.^3,4 More recently, the effects on stability and structural aspects of a protein of the presence of ionizable pair buried in the hydrophobic core of a native protein are investigated.^5–14 Though, the interacting buried ionizable pairs (an acidic and basic residue) are not omnipresent in nature, they are important in many biological process like catalysis, electron transfer and proton transport etc. García-Moreno and coworkers have artificially engineered a Glu/Lys (V23E/L36K) pair in the stable hydrophobic core of a hyperstable variant of staphylococcal nuclease (Δ+PHS) protein (PDB ID: 3BDC).²⁰

Their work is aimed to understand the thermodynamic and structural consequences in the artificially engineered protein. They try to address key issues, namely, (i) whether the favorable coulomb interaction between the polarizable charge pair can compensate for the unfavorable self-energy experienced by the two buried charged acidic and basic residues in the dry core of the engineered protein and (ii) if the scaffold protein is stable in a wide range of environmental conditions like temperature, pH etc., then it can be used to develop a novel class of catalysts, without the need for specialized structural modifications. For this system, it is observed that the global structure of the protein is unaffected by the introduction of the internal buried polarizable pair, however, alterations to the microenvironment around the ion pair are observed with Gly-20 and Thr-62 rearranging themselves to interact with Lys-36 and Glu-23, respectively (see Fig. 1).

Fig. 1 Microenvironment of the V23E/L36K variant of Δ+PHS protein showing several important core residues, Gly-20, Asp-21, Glu-23, Arg-35, Lys-36 and Thr-62. The structure shown here is the representative structure of the dominant cluster based on RMS cluster analysis from 25 ns CpHMD simulation conducted at pH 8.

The pK_a of a particular residue control the protonation/deprotonation state of that particular titratable amino acid residue at a particular pH of the environment. The pK_a of a particular residue is also influenced by the interaction of other neighboring residues. Conventional molecular dynamics (cMD) simulations are not equipped to account for the interplay between conformation and its protonation state (protonation/deprotonation). Various theoretical techniques have been developed to consider pH explicitly in MD simulation and successfully applied in various types' proteins.^21–43 The constant pH MD simulation (CpHMD) method employed in the present work was developed by Mongan et al., which uses an implicit solvent model, namely, generalized Born (GB).²⁵

The aim of the present manuscript is to analyze and validate whether the presently adopted CpHMD calculation with GB implicit model for the solvent water can reproduce the crystallographic conformation and can capture the minor rearrangements of Gly-20 and Thr-62 residues in the microenvironment of the doubly mutated (V23E/L36K) variant, when one starts from the parent form, Δ+PHS (PDB ID: 3BDC). It will be also interesting to study the hydrogen bond dynamics in the microenvironment by constructing a suitable correlation function. Present study will also help to validate and augment recent experimental findings.²⁰

Theoretical methodology

Constant pH molecular dynamics simulation

In the present study, we have considered a mutated variant of staphylococcal nuclease (SNase), Δ+PHS (PDB ID: 3BDC).⁹ Relative to the wild type SNase, Δ+PHS is hyper-stable and acid-resistant variant containing five substitutions (G50F, V51N, P117G, H124L and S128A).^44,45 To generate the variant containing a charged pair studied by Garcia-Moreno and coworkers, a Glu/Lys pair (V23E/L36K) is then virtually introduced into the hydrophobic core of the Δ+PHS form by mutation using the Maestro software package.⁴⁶ The doubly mutated structure is then used for CpHMD simulation. The CpHMD method employs MD with a generalized Born (GB) implicit model for the solvent.²⁵ During the course of the simulation, the MD simulation is periodically halted, and a Monte-Carlo (MC) step is taken. In the MC step a titratable residue (Glu, Asp, His, Lys, Tyr and Cys) is considered randomly to change its protonation state. The pK_a is calculated with respect to a reference compound for the residue of interest. The reference compounds are the isolated titratable residues solvated in water and the reference pK_a values are 3.8 for Asp, 4.3 for Glu, 6.8 for His, 9.6 for Tyr and 10.5 for Lys. The pK_a for a desired residue is calculated by using the transition energy corresponding to the MC step, expressed as,

ΔG = k_BT(pH − pK_a,ref)ln 10 + ΔG_elec − ΔG_elec,ref (1)

where, T is absolute temperature, pH is the solvent pH, k_B is the Boltzmann constant, pK_a,ref is the pK_a of the reference compound, ΔG_elec and ΔG_elec,ref are the electrostatic energy change for protonation state change of the desired residue and the reference compound, respectively. The transition energy for the pK_a calculation is based on the GB electrostatics applied in the MD simulation. The titratable residues are selected randomly for titration and acceptance of the change in protonation state is solely governed by the Metropolis criterion. If the MC move is accepted based on the Metropolis criterion, then the MD run is continued with the changed protonation state for the titratable residue. On the other hand if the MC move is rejected based on the Metropolis criterion, then the MD run is continued with the present protonation state for that residue.

The CpHMD method with MC steps for titration is implemented in the AMBER12 suite of program through the Sander molecular dynamics engine.⁴⁷ All the MD simulations are carried out using the AMBER99SB force field.⁴⁸ A GB implicit solvent model (igb = 5)^49–51 with salt concentration of 0.1 M is employed using a 30 Å cutoff value for nonbonded interactions and for calculation of effective Born radii. The SHAKE algorithm is used to constrain all the bonds involving hydrogen.⁵² The Berendsen thermostat with a relaxation time of 2 ps is used to maintain the temperature at 300 K.⁵³ MD time step is taken as 2 fs and a period of 10 fs of MD separates the MC titrations. The system is heated for 1 ns in NPT ensemble and then passes through a 2 ns NVT equilibration run before the production run of 10 ns over a pH range of 3 to 11 at increments of 0.5 pH unit. Simulations at extreme acidic pH (<3) and basic pH (>11) are not considered in the present study, as the doubly substituted variant shows denaturation in at extreme pH.²⁰ All visualizations are carried out by employing the VMD program.⁵⁴

Quantum mechanical calculation

Full geometry optimizations are carried out for the glutamate–lysine complex using B3LYP/6-311++G(d,p) level of theory. B3LYP hybrid functional consists of Becke's three parameters non-local exchange (B3) and Lee–Yang–Parr (LYP) non-local correlation.⁵⁵ A pseudo Newton–Raphson based algorithm is used to carry out geometry optimization in the multidimensional potential energy surface for the glutamate–lysine complex with various initial structures. The major drawback of this geometry optimization methodology is to guess good initial structure, which might converge to a minimum energy structure during full geometry search. For this a large number of (ten) initial conformations are considered in the present most stable structure search. Conductor like screening model is applied to model the solvent continuum.⁵⁶ All quantum mechanical calculations are carried out by employing GAMESS suite of ab initio program system.⁵⁷

Results and discussion

Conformation and geometry

As an initial test of performance CpHMD simulations relative to experiment, various distances between the side chains of Glu-23 and Lys-36 are calculated over the 10 ns trajectories and are compared to the crystallographic structure for the V23E/L36K variant [PDB ID: 3NHH]. The various average calculated distances between Nζ of Lys-36 and Oε1 and Oε2 of Glu-23 and between Cε of Lys-36 and Cδ of Glu-23 are shown in Table 1. Root-mean-square (RMS) clustering analysis are performed on the trajectories obtained at each pH level to aid in visualizing the dominant conformations sampled by the systems. The Glu-23 and Lys-36 residues are also shown in a representative structure for the dominant cluster obtained from cluster analysis of trajectories at pH 8 in Fig. 1, by taking into consideration that the experimental crystal structure for the V23E/L36K variant (PDB ID: 3NHH) is reported at a pH of 8. The various geometrical parameters are also calculated from the optimized geometry of glutamate–lysine complex at B3LYP/6-311++G(d,p) level and are shown in Table 1. Various geometrical parameters are also calculated for protonated and deprotonated states by employing conventional MD (cMD) simulation with the same simulation parameters from 25 ns production run and the values are also shown in Table 1. It is interesting to observe from Table 1 that various distances between the buried charged Glu/Lys pair at pH 8 based on CpHMD simulation and quantum mechanical calculation are very close to that of the experimentally reported crystal structure taken at pH 8, demonstrating that the CpHMD simulation can reliably reproduce the various critical distances between the buried ionizable charged pair.

Table 1 Various average distances (Å) between the Glu-23/Lys-36 pair at different level of pH for V23E/L36K variant of Δ+PHS averaged over 10 ns CpHMD trajectories. Same distances are shown for the crystallographic structure (PDB ID: 3NHH), which is recorded at pH 8

pH	Oε1–Nζ	Oε2–Nζ	Cδ–Cε
a Values in the parentheses are for 25 ns run.b Values are based on conventional MD.
3	6.12	6.35	5.90
3.5	5.56	5.82	5.88
4	4.75	5.12	5.20
4.5	5.20	5.35	5.50
5	4.65	4.81	5.05
5.5	4.45	4.61	4.87
6	4.64	4.76	5.15
6.5	3.76	4.05	4.71
7	6.35	4.90	5.37
7.5	3.56	3.13	4.30
8^a	2.98 (2.99)	3.14 (3.16)	4.37 (4.32)
8.5	2.85	3.75	4.53
9.0	2.85	4.08	4.81
9.5	3.10	3.22	4.75
10	3.35	3.52	4.80
10.5	4.61	4.77	5.05
11	5.25	5.53	5.76
Protonated^b	4.67	4.87	5.22
Deprotonated^b	4.42	4.69	5.06
Quantum mechanical (6-311++G(d,p))	2.59	3.02	4.13
Crystal structure	2.56	3.13	4.04

---

RMSD and RMSF in the microenvironment and surface

To assess the stability of CpHMD simulations we have extended the simulations up to 25 ns at pH 8, the pH value for which the double mutant crystal structure is reported. The root mean square displacement (RMSD) of the backbone of the full scaffold protein is calculated from the CpHMD simulation with respect to the crystal structure of V23E/L36K variant (PDB ID: 3NHH) and the results are shown in Fig. 2 respectively. The whole backbone fluctuation is shown marked as ‘a’ (black line) for both plots. We have also considered the microenvironment within 10 Å and 5 Å of residue 36 (Lys-36) for RMSD calculations to assess convergence in the simulation. The backbone RMSD of the protein residues within 10 Å and 5 Å of Lys-36 are shown marked as ‘b’ (red line) and ‘c’ (blue line), respectively, in the figure. It is clear from the figure that the simulation is converged. Superposition of the simulated structure with the crystal structures of V23E/L36K variant (PDB ID: 3NHH) is shown in Fig. 3 for comparison.

Fig. 2 RMSD (Å) with respect to the crystal structure of V23E/L36K variant (PDB ID: 3NHH) for the backbone from a CpHMD simulation at pH 8. The black line (marked as ‘a’) in each figure represents the RMSD for the entire protein. The red line (marked as ‘b’) and blue line (marked as ‘c’) represent the RMSD fluctuations for residues residing within 10 and 5 Å of Lys-36.

Fig. 3 Superposition of the simulated structure (cyan) for V23E/L23K variant at pH 8 using CpHMD with the crystal structure of V23E/L23K (red, PDB ID: 3NHH) after a 25 ns simulation. The simulated structure considered here is the representative structure of the dominant cluster based on RMS cluster analysis.

We have also calculated backbone root mean square fluctuations (RMSF) over the trajectory for all the residues, as shown in Fig. 4. It is worthwhile to mention that the RMSF is calculated on trajectory that is aligned to the all backbone atoms. RMSF of the residues within 10 and 5 Å of Lys-36 are also shown in the figure. It is clear from the RMSF calculation, that the residues within 5 Å of Lys-36 fluctuate least (<2 Å), with the exception of Thr-62 (RMSF = 2.2 Å). The high fluctuations of Thr-62 can be attributed to a rearrangement of side chain (see details in the next section). The other residues which displayed more than 2 Å RMSF are residue number 1–2 (terminal residues), 27–30, 69–71, 77–87, 95–96, 114–141. All of these residues are beyond the region within 5 Å of Lys-36, with a majority of them residing at the surface. Based on RMSD and RMSF calculation it is clearly observed that the core of the scaffold protein does not change upon introducing a charge pair deep inside the hydrophobic core (see also Fig. 3).

Fig. 4 The RMSF fluctuations (black solid line) in Å for the backbone protein of the V23E/L36K variant from a 25 ns CpHMD simulation at pH 8. The open red circles and blue squares represent the residues within 10 and 5 Å of Lys-36, respectively.

Hydrogen bonds in the microenvironment

We have measured the hydrogen bond formation of various residues within 10 Å around residue 36 (Lys-36) over the 25 ns trajectory based on CpHMD simulation at pH 8, and compared this analysis with the report of hydrogen bonds available experimentally for the Δ+PHS system (see Table 2). The distance and angle cut off for hydrogen bond calculations are taken as 3 Å (0.3 nm) and 30°, respectively. All the hydrogen bonds observed experimentally for the Δ+PHS system are preserved in our simulations of the doubly mutated variant, with the exception between the K70-D95 pair. These observations clearly suggest that the global conformation does not change significantly by introducing a charge pair in the hydrophobic core.

Table 2 Hydrogen bonding in the micro-environment (within 10 Å from the residue K36) of the V23E/L36K variant of Δ+PHS from a 25 ns CpHMD simulation conducted at pH 8

Amino acids pair	Δ+PHS^a	V23E/L36K variant ^b
a Values are taken from NMR experiment obtained by the Garcia-Moreno lab (PDB ID: 3BDC).^9b Values in the parentheses are the percentage (%) of time a particular hydrogen bond survives in the 25 ns run.
T41–D21	Yes	Yes (65)
E52–E43	Yes	Yes (38)
E57–Y54	Yes	Yes (1)
D19–D21	Yes	Yes (1)
K70–D95	Yes	No
K71–D95	Yes	Yes (2)
E57–E57	Yes	Yes (<1)
D83–Y85	Yes	Yes (10)
D83–R87	Yes	Yes (68)
Y85–R87	Yes	Yes (<1)
D21–R35	Not explored in experiment	Yes (16)
E23–K36	—	Yes (19)
G20–K36	—	Yes (20)
E23–T62	—	Yes (10)

---

Significant changes instilling a hydrogen bond network in the vicinity of the mutated E23/K36 pair is observed in our present studies. In the CpHMD simulations, a minor rearrangement of the backbone carbonyl oxygen of Gly-20 and the side chain hydroxyl group of Thr-62 occur that allow these moieties to interact with the Lys-36 and Glu-23, respectively, in V23E/L36K. Rearrangement of the backbone carbonyl oxygen of G20 and the side chain hydroxyl group of T62 in the micro-environment of the E23 and K36 charged pair are shown in Fig. 5. The critical distances for these interactions are also provided in the Table 3 for our simulation, showing good agreement with the experimental crystal structure. Form these interactions, three extra hydrogen bonds are observed for the double mutant that are not present in Δ+PHS (see Table 2), namely, between E23–K36, E23–T62 and G20–K36 pairs in the V23E/L36K variant. All these hydrogen bonds are solely attributed due to the introduction of a Glu-23/Lys-36 pair in the hydrophobic core.

Fig. 5 Rearrangement of the backbone carbonyl oxygen of G20 and the side chain hydroxyl group of T62 in the micro-environment of the E23 and K36 charged pair seen in a 25 ns CpHMD simulation conducted at pH 8. The residues are also shown. In the inset, the residues of the crystal structures (licorice) and the representative structure of the dominant cluster from CpHMD simulation (ball and stick) show good agreement.

Table 3 Various distances (Å) between Gly-20 and Lys-36 as well as between Thr-62 with Glu-23 for the V23E/L36K variant over a 25 ns CpHMD simulation conducted at pH 8

Structure	O (Gly-20)–Nζ (Lys-36)	Hydroxyl O (Thr-62)–Oε2 (Glu-23)
CpHMD at pH 8	2.97	4.14
Crystal structure	2.93	3.97

---

Time dependent behavior of hydrogen bond is also important.^58,59 The hydrogen bond auto correlation function, C_HB(t) is also calculated as

where, n is the total number of hydrogen bond forms in the microenvironment, h(0) and h(t) are the hydrogen bond at time 0 and t, respectively. The decay of auto correlation function in the microenvironment (within 10 Å of K36) is shown in Fig. 6. The decay of the hydrogen bond auto correlation function is excellently fitted (χ² = 10⁻⁴ and R = 0.98) with single exponential curve having time constant of 2282 ps. The best fitted expression is

Fig. 6 Time dependence of the hydrogen bond correlation function, C_HB(t) for all the hydrogen bond formed in the micro-environment (within 10 Å from the residue K36) of the V23E–L36K variant.

We observe a very long hydrogen bond correlation. The hydrogen bond lifetime is also calculated to be 716 ps. These observations clearly state the stability of the core of the scaffold protein.

Salt bridge interactions

The distance between the interacting charged groups is very important for salt bridge formation, and usually it should be around 4 Å to be called as salt bridge interaction.⁶⁰ Both electrostatic interactions and hydrogen bonding between the charged groups are responsible for the formation of salt bridge. Salt bridge interactions not only play a vital role in the thermodynamic stability or instability of a particular folded or unfolded state of a protein but are also important in fine tuning the structure of a protein for its optimum function.^61–64 Present CpHMD simulations provide strong evidence of salt-bridge formation between Asp-21/Arg-35, Glu-23/Lys-36, Arg-35/Tyr-113 and Asp-40/Lys-110 pairs. Based on the CpHMD calculation, the average distance between the interacting groups of these four pairs is less than 4 Å. It is interesting to observe that, two surface residues, Lys-110 and Tyr-113 are also involved in a salt bridge interactions. All these strong salt bridge interactions lead to stability of the core, as well as of the whole scaffold protein. The strong salt bridge interaction between the mutated pair, E23–K36, coupled with hydrogen bonding between E23–T62 and G20–K36 pairs (see Fig. 5) in the V23E/L36K variant, may help to compensate for the unfavorable self-energy experienced by these residues, when buried in the hydrophobic core.

Conclusions

The adopted CpHMD calculation with a GB implicit model for the solvent water is able to reproduce the crystallographic conformation of the doubly mutated (V23E/L36K) variant, when one starts from the parent form, Δ+PHS. It is observed that the average geometrical parameters based on CpHMD calculation at pH 8 nicely match to those reported for the crystal structure of the V23E/L36K variant. The geometrical parameters calculated at B3LYP/6-311++G(d,p) level are also very good in agreement with that of the crystal structure. Strong hydrogen bond formation and salt bridge interaction in the microenvironment around the Glu-23/Lys-36 pair stabilize the core of the V23E/L36K variant. We do observe very long hydrogen bond correlation and long hydrogen bond life time (718 ps) in the microenvironment. We also observe very good single exponential fitting of the decay of the hydrogen bond auto correlation function with a time constant of 2282 ps. Small rearrangements of the backbone carbonyl oxygen of Gly-20 and the side chain hydroxyl group of Thr-62 are also observed, and these changes allow them to interact with Lys-36 and Glu-23, respectively in V23E/L36K. This is also consistent with the experiment. We have observed three new hydrogen bonds in comparison to the Δ+PHS form (see Table 2), namely between E23–K36, E23–T62 and G20–K36 pairs in the V23E/L36K variant. Based on the present study, no significant change in the conformation of V23E/L36K with respect to the parent form Δ+PHS is noticed. Overall, the CpHMD calculation with implicit solvent model can efficiently capture fine rearrangements of residues in a protein with a charged pair buried in its hydrophobic core and the method can be very useful for blind prediction study.

Acknowledgements

This work was supported in part by the IUSSTF and UCSD. The Author would like to acknowledge Prof. Bertrand E. Garcia-Moreno for sharing his results on the V23E/L36K variant (ref. 20) and Prof. J. A. McCammon and Dr Ptarick Blachly of UCSD for helpful discussions. Dr(s) B. N. Jagatap, T. Bandyopadhyay, A. K. Samanta are acknowledged for their encouragement. Author declares no conflict of interest.

References

1. C. Pace, B. Shirley, M. McNutt and K. Gajiwala, FASEB J., 1996, 10, 75–83.
2. G. Rose, P. Fleming, J. Banavar and A. Maritan, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 16623–16633.
3. M. S. Chimenti, et al., Structure, 2012, 20, 1071–1085.
4. C. A. Castañeda, et al., Proteins: Struct., Funct., Bioinf., 2009, 77, 570–588.
5. C. von Ballmoos, A. Wiedenmann and P. Dimroth, Annu. Rev. Biochem., 2009, 78, 649–672.
6. R. H. Fillingame, C. M. Angevine and O. Y. Dmitriev, FEBS Lett., 2003, 555, 29–34.
7. A. Aksimentiev, I. A. Balabin, R. H. Fillingame and K. Schulten, Biophys. J., 2004, 86, 1332–1344.
8. S. Iwata, et al., Science, 1998, 281, 64–71.
9. T. Tsukihara, et al., Science, 1996, 272, 1136–1144.
10. T. Tsukihara, et al., Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 15304–15309.
11. T. K. Harris and G. J. Turner, IUBMB Life, 2002, 53, 85.
12. A. Warshel, Proc. Natl. Acad. Sci. U. S. A., 1978, 75, 5250–5254.
13. W. R. Cannon and S. J. Benkovic, J. Biol. Chem., 1998, 273, 26257–26260.
14. A. K. Pathak, Biopolymers, 2015, 103, 148.
15. S. Deechongkit, H. Nguyen, P. E. Dawson, M. Gruebele and J. W. Kelly, Nature, 2004, 403(45), 101–105.
16. S. T. Bayley, Trans. Faraday Soc., 1951, 47, 509.
17. J. J. Dwyer, et al., Biophys. J., 2000, 79, 1610–1620.
18. D. G. Isom, B. R. Cannon, C. A. Castañeda, A. Robinson and E. B. García-Moreno, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 17784–17788.
19. A. K. Pathak, Mol. Phys. 2015, 113, 3345–3352.
20. A. C. Robinson, C. A. Castañeda, J. L. Schlessman and B. García-Moreno, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 11685–11690.
21. Y. Y. Sham, Z. T. Chu and A. Warshel, J. Phys. Chem. B, 1997, 101, 4458–4472.
22. K. M. Merz, J. Am. Chem. Soc., 1991, 113, 3572–3575.
23. J. H. Jensen, H. Li, A. D. Robertson and P. A. Molina, J. Phys. Chem. A, 2005, 109, 6634–6643.
24. D. Riccardi, P. Schaefer, Y. Yang, H. B. Yu, N. Ghosh, X. Prat-Resina, P. Konig, G. H. Li, D. G. Xu, H. Guo, M. Elstner and Q. Cui, J. Phys. Chem. B, 2006, 110, 6458–6469.
25. J. Mongan, D. A. Case and J. A. J. McCammon, J. Comput. Chem., 2004, 25, 2038–2048.
26. M. O. Kim and J. A. McCammon, Biopolymers,  DOI:10.1002/bip.22702.
27. M. Machuqueiro and A. M. Baptista, Proteins: Struct., Funct., Bioinf., 2008, 72, 289–298.
28. M. Dlugosz and J. M. Antosiewicz, Chem. Phys., 2004, 302, 161–170.
29. M. Dlugosz, J. M. Antosiewicz and A. D. Robertson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 021915.
30. A. M. Baptista, P. J. Martel and S. B. Petersen, Proteins: Struct., Funct., Genet., 1997, 27, 523–544.
31. U. Borjesson and P. H. Hunenberger, J. Chem. Phys., 2001, 114, 9706–9719.
32. U. Borjesson and P. H. J. Hunenberger, J. Phys. Chem. B, 2004, 108, 13551–13559.
33. M. S. Lee, F. R. Salsbury and C. L. Brooks, Proteins: Struct., Funct., Bioinf., 2004, 56, 738–752.
34. J. Khandogin and C. L. Brooks, Biophys. J., 2005, 89, 141–157.
35. J. Khandogin, J. H. Chen and C. L. Brooks, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 18546–18550.
36. S. L. Williams, C. A. F. de Oliveira and J. A. McCammon, J. Chem. Theory Comput., 2010, 6, 560–568.
37. Y. L. Meng and A. E. Roitberg, J. Chem. Theory Comput., 2010, 6, 1401–1412.
38. E. Alexov, E. L. Mehler, N. Baker, A. M. Baptista, Y. Huang, F. Milletti, J. E. Nielsen, D. Farrell, T. Carstensen, M. H. M. Olsson, J. K. Shen, J. Warwicker, S. Williams and J. M. Word, Proteins: Struct., Funct., Bioinf., 2011, 79, 3260–3275.
39. E. J. Arthur, J. D. Yesselman and C. L. Brooks III, Proteins: Struct., Funct., Bioinf., 2011, 79, 3276–3286.
40. M. H. M. Olsson, Proteins: Struct., Funct., Bioinf., 2011, 79, 3333–3345.
41. J. A. Wallace, Y. Wang, C. Shi, K. J. Pastoor, B.-L. Nguyen, K. Xia and J. K. Shen, Proteins: Struct., Funct., Bioinf., 2011, 79, 3364–3374.
42. S. G. Itoh, A. Damjanović and B. R. Brooks, Proteins: Struct., Funct., Bioinf., 2011, 79, 3420–3436.
43. S. L. Williams, P. G. Blachly and J. A. McCammon, Proteins: Struct., Funct., Bioinf., 2011, 79, 3381–3388.
44. D. G. Isom, B. R. Cannon, C. A. Castaneda, A. Robinson and B. E. Garcia-Moreno, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 17784–17788.
45. D. Bashford, D. A. Case, C. Dalvit, L. Tennant and P. E. Wright, Biochemistry, 1993, 32, 8045–8056.
46. Maestro, version 9.3, Schrödinger, LLC, New York, NY, 2012.
47. D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz, B. Roberts, S. Hayik, A. Roitberg, G. Seabra, J. Swails, A. W. Goetz, I. Kolossváry, K. F. Wong, F. Paesani, J. Vanicek, R. M. Wolf, J. Liu, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M.-J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, R. Salomon-Ferrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko and P. A. Kollman, AMBER 12, University of California, San Francisco, 2012.
48. J. Chen, W. Im and C. L. Brooks III, J. Am. Chem. Soc., 2006, 128, 3728–3736.
49. A. Onufriev, D. A. Case and D. Bashford, J. Comput. Chem., 2002, 23, 1297–1304.
50. A. Onufriev, D. Bashford and D. A. Case, J. Phys. Chem. B, 2000, 104, 3712–3720.
51. A. Onufriev, D. Bashford and D. A. Case, Proteins: Struct., Funct., Bioinf., 2004, 55, 383–394.
52. J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
53. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
54. W. Humphrey, A. Dalke and K. Schulten, VMD – Visual Molecular Dynamics, J. Mol. Graphics, 1996, 14(1), 33–38.
55. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
56. K. Baldridge and A. Klamt, J. Chem. Phys., 1997, 106, 6622–6633.
57. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis and J. A. J. Montgomery, J. Comput. Chem., 1993, 14, 1347.
58. A. Luzar and D. Chandler, Nature, 1996, 379, 55.
59. A. Luzar and D. Chandler, Phys. Rev. Lett., 1996, 76, 928.
60. S. Kumar and R. Nussinov, ChemBioChem, 2002, 3, 604–617.
61. Z. S. Hendsch, T. Jonsson, R. T. Sauer and B. Tidor, Biochemistry, 1996, 35, 7621–7625.
62. R. Norkel, F. Sheinerman, D. Petrey and B. Honig, Protein Sci., 2001, 10, 2147–2161.
63. Z. S. Hendsch and B. Tidor, Protein Sci., 1994, 3, 211–226.
64. F. B. Sheinerman and B. Honig, J. Mol. Biol., 2002, 318, 161–177.

---