Flexibility of Catalytic Zinc Coordination in Thermolysin and HDAC8: A

Dec 2, 2009 - Tamames , B.; Sousa , S. F.; Tamames , J.; Fernandes , P. A.; Ramos , M. J. ... Sousa , S. F.; Fernandes , P. A.; Ramos , M. J. J. Am. C...
0 downloads 0 Views 2MB Size
J. Chem. Theory Comput. 2010, 6, 337–343

337

Flexibility of Catalytic Zinc Coordination in Thermolysin and HDAC8: A Born-Oppenheimer ab Initio QM/MM Molecular Dynamics Study Ruibo Wu,†,‡ Po Hu,† Shenglong Wang,† Zexing Cao,‡ and Yingkai Zhang*,† Department of Chemistry, New York UniVersity, 100 Washington Square East, New York 10003, and Department of Chemistry and State Key Laboratory of Physical Chemistry of Solid Surfaces, College of Chemistry and Chemical Engineering, Xiamen UniVersity, Xiamen 361005, China Received October 7, 2009

Abstract: The different coordination modes and fast ligand exchange of zinc coordination has been suggested to be one key catalytic feature of the zinc ion which makes it an invaluable metal in biological catalysis. However, partly because of the well-known difficulties for zinc to be characterized by spectroscopy methods, evidence for dynamic nature of the catalytic zinc coordination has so far mainly been indirect. In this work, Born-Oppenheimer ab initio Quantum Mechanical/Molecular Mechanical (QM/MM) molecular dynamics (MD) simulation has been employed, which allows for a first-principle description of the dynamics of the metal active site while properly including effects of the heterogeneous and fluctuating protein environment. Our simulations have provided direct evidence regarding inherent flexibility of the catalytic zinc coordination shell in thermolysin (TLN) and histone deacetylase 8 (HDAC8). We have observed different coordination modes and fast ligand exchange during the picosecond’s time scale. For TLN, the coordination of the carboxylate group of Glu166 to zinc is found to continuously change between monodentate and bidentate manner dynamically, while for HDAC8, the flexibility mainly comes from the coordination to a nonamino acid ligand. Such distinct dynamics in the zinc coordination shell between two enzymes suggests that the catalytic role of zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of the amide bond. Meanwhile, considering that such Born-Oppenheimer ab initio QM/MM MD simulations are very much desired but are widely considered to be too computationally expensive to be feasible, our current study demonstrates the viability and powerfulness of this state-of-the-art approach in simulating metalloenzymes.

Zinc is relatively abundant in biological materials. Approximately 10% of the total human proteome have been identified to bind with zinc in vivo from a bioinformatics investigation,1 and they play very crucial roles in all forms of life.2-6 For mononuclear zinc enzymes, a typical metal coordination environment contains three amino acid side chain ligands (His, Glu, Asp, and Cys) and one or two small * To whom correspondence should be addressed. E-mail: [email protected]. Telephone: (212) 998 7882. † New York University. ‡ Xiamen University.

molecules.3,7,8 The flexibility of zinc coordination, which allows different coordination modes and fast ligand exchange, has been suggested to be one key catalytic feature of the zinc ion which makes it an invaluable metal in biological catalysis.9 However, partly because of the well-known difficulties for zinc to be characterized by spectroscopy methods,10,11 evidence for dynamic nature of the catalytic zinc coordination has so far mainly been indirect coming from determined X-ray crystal structures5,12 and geometry optimizations of model complexes with electronic structure methods.13,14

10.1021/ct9005322  2010 American Chemical Society Published on Web 12/02/2009

338

J. Chem. Theory Comput., Vol. 6, No. 1, 2010

A classical example of the flexible zinc coordination is thermolysin (TLN), one of the best experimentally characterized zinc proteases with one Glu and two His amino acid side chains as ligands.15 Extensive structural studies have indicated that besides the simple tetrahedron coordination, its catalytic zinc ion can also be five- or six-coordinated16-19 especially for the carboxylate of Glu/Asp whose coordination is highly flexible. The different coordination mode of carboxylate to zinc with a continuous range between monodentate and bidentate manner has been observed in model complexes as well as in other zinc enzymes with one Glu or Asp as ligands.12 It has been suggested that the coordination of carboxylate to zinc could be dynamic in the catalytic process, known as carboxylate shift, which has been suggested to be important in the function of zinc enzymes.13,14,20-22 Intriguingly, for zinc-dependent histone deacetylases (HDACs), in which the catalytic zinc ion is coordinated to one His and two Asp residues, only the monodentate mode has been observed for the carboxylate coordination in several crystal structures.23-26 HDACs, which catalyze the removal of acetyl group from histone tails, have emerged to be critical in gene regulation and are among the most promising targets for the development of antitumor therapeutics. For example, the recently FDA approved anticancer drug SAHA is an HDAC inhibitor which is directly coordinated to the zinc.27,28 To provide deep insights into the dynamics and flexibility of the zinc catalytic site, which would be essential in characterizing their catalytic mechanisms and rational design of novel inhibitors for zinc enzymes, we have carried out density functional theory (DFT) QM/MM Born-Oppenheimer molecular dynamics (BOMD) simulations on TLN and HDAC8. Although semiempirical QM/MM BOMD simulations of some zinc-dependent enzymes have been carried out,29-32 one main concern is the accuracy and transferability of the semiempirical QM Hamiltonian in describing the zinc coordination shell and the enzyme reaction profile. Often, some parameters of the semiempirical Hamiltonian need to be reoptimized for specific systems and reactions to obtain reasonable results,30,32-35 which significantly reduce its transferability and predictability. Thus, in our current DFT QM/MM molecular dynamics simulations, the zinc ion and its ligands are treated by B3LYP hybrid functional with a Stuttgart ECP/basis set36 for the zinc atom and a 6-31G* basis set for all other QM atoms. Our employed theoretical approach, which has recently been demonstrated to be feasible and powerful in several enzymatic studies,37-40 allows for a first-principle description of the dynamics of the metal active site while properly including effects of the heterogeneous and fluctuating protein environment.

Computational Methods I. Preparation of Simulation Systems. As illustrated in Figure 1 and listed in Tables 1 and 2, four TLN and two HDAC8 simulation systems have been prepared on the basis of their respective experiment crystal structures: TLNa, a substrate-free structure of TLN (pdb code: 1LNF16); TLNb, a TLN-inhibitor complex which mimics the reactant state (pdb code: 1ZDP17); TLNc, a TLN-inhibitor complex which

Wu et al.

Figure 1. Illustration of various zinc coordination shells in our simulated TLN and HDAC8 enzyme complexes.

mimics the transition state (pdb code: 2TMN18); TLNd, a TLN-product complex (pdb code: 3TMN19); model 1, a complex of the Y306F mutant of HDAC8 and its substrate (pdb code: 2 V5W25); model 2, an HDAC8-SAHA complex (pdb code: 1T6924). Representative active site structures of TLN and HDAC8 enzymes are illustrated in Figure 2.

Thermolysin and HDAC8 Catalytic Zinc Coordination

J. Chem. Theory Comput., Vol. 6, No. 1, 2010 339

Table 1. Interaction Distances in the Zinc Coordination Sphere of Thermolysina CNb

Zn-N (H142) d1

Zn-N (H146) d2

Zn-Oe1(E166) d3

Zn-Oe2(E166) d4

d5

d6

6 5.6 ( 0.4 4 4.0 ( 0.1 4 5.3 ( 0.4 4 4.5 ( 0.4

1.980 2.063 ( 0.063 2.141 2.031 ( 0.051 2.119 2.043 ( 0.051 2.119 1.991 ( 0.050

1.987 2.152 ( 0.076 2.096 2.066 ( 0.060 2.102 2.235 ( 0.116 2.117 2.078 ( 0.062

2.378 2.216 ( 0.289 2.002 2.053 ( 0.062 2.148 2.940 ( 0.409 2.167 2.173 ( 0.272

2.242 2.150 ( 0.069 2.955 3.117 ( 0.140 2.887 2.131 ( 0.335 2.800 2.466 ( 0.415

2.384g 2.186 ( 0.093g 2.414h 2.330 ( 0.060h 2.792i 2.220 ( 0.223i 2.129j 2.089 ( 0.071j

2.276g 2.158 ( 0.086g n/a n/a 2.065i 2.083 ( 0.128i n/a n/a

models c

expt (1LNF) TLNa expt (1ZDP)d TLNb expt (2TMN)e TLNc expt (3TMN)f TLNd

a The flexible distances are highlighted in bold (n/a, not applicable). Å are used for Zn-O, 2.40 and 2.65 for Zn-S, respectively. These values are chosen on the basis of very recent analysis on databases of zinc enzyme structures database.5 b CN means the coordination number. CN is 1 if Zn-N e 2.15 Å, equals 0 if Zn-N g 2.40 Å, and is a linear scalar between 0 and 1 if Zn-N is between 2.15 and 2.40 Å. Similarly, the values of 2.20 and 2.60. c Ref 15. Resolution: 1.7 Å. d Ref 16. Resolution: 1.7 Å. e Ref 17. Resolution: 1.6 Å. f See from ref 18. Resolution: 1.7 Å. g It is Zn-O distance between Zn2+ and crystal water. h It is Zn-S distance between Zn2+ and inhibitor. i It is Zn-O distance between Zn2+ and phosphate group in inhibitor. j It is Zn-O distance between Zn2+ and crystal water.

Table 2. Interaction Distances in the Zinc Coordination Sphere of HDAC8 models

CN

Zn-O (D178) r1

Zn-N (D180) r2

Zn-O (D267) r3

Zn-(water/SAHA) r4

Zn-(sub./SAHA) r5

expt (2 V5W)a model 1 (Y306F) expt (1T69)b model 2 (SAHA)

5 4.4 ( 0.4 5 4.8 ( 0.2

1.955 1.938 ( 0.047 2.131 1.992 ( 0.056

2.074 1.996 ( 0.049 1.789 2.154 ( 0.078

1.947 2.002 ( 0.085 1.908 1.969 ( 0.052

2.238 2.422 ( 0.219 1.953 2.247 ( 0.111

1.994 2.103 ( 0.064 1.975 2.171 ( 0.083

a

Reference 25. Resolution: 2.0 Å.

b

Reference 24. Resolution: 2.9 Å.

For each simulation system, its initial structure was prepared on the basis of its respective crystal structure with molecular modeling and dynamics simulations. First, missing residues were added with the Swiss-PdbViewer.41 Protonation states of charged residues in the enzyme complex were determined via H++ program42 and by carefully examining local hydrogen bond networks. In all TLN models, His142 and His146 were identified as singly protonated, whereas His231 was doubly protonated, which are the same as previous studies.43,44 For HDAC8 models, His180 was determined as singly protonated on δ site, while His142 and His143 were determined as singly protonated on ε site. Then, each prepared system was solvated into a rectangular box with a 10 Å buffer distance between the solvent box wall and the nearest solute atom. Finally, to neutralize each simulation system, one to six sodium ions were added at the protein surface by employing the Amber Tool. All sodium ions were located more than 17 Å away from the zinc active site. The resulting simulation system was about 45 000 atoms for each model. The Zn2+ ion was modeled with the Stote’s scheme.45 Considering that the zinc coordination shell is very difficult to be well described by a molecular mechanical force field,46-49 ∼1500 kcal/mol harmonic restraint was placed on ∼15 atoms in the zinc bind site to retain the zinc coordination structure during the MM equilibrations. The rest of protein and solvent molecules were first minimized, and then more than 3 ns MD simulations were carried out for each system with periodic boundary condition. A time step of 2 fs was used. Berendsen thermostat method50 was used to control the system temperature at 300 K. All MD simulations were performed with AMBER10 molecular dynamics package.51 The Amber99SB52-54 force field for the protein and TIP3P model55 for water molecules were employed, and the force field parameters for substrate in HDAC8 and inhibitors in these simulation models were generated from AMBER GAFF force field56 via AMBER tools. The SHAKE algorithm57 was applied to constrain all

bonds involving hydrogen atoms with tolerance of 10-5 Å, and 12 Å cutoffs were used for both van der Waals and electronic interactions. II. Born-Oppenheimer ab Initio QM/MM Molecular Dynamic Simulations. Considering that the trajectory was very stable after 2 ns for all these classical MD simulations, the resulting snapshot after 3 ns MD simulation was employed as the initial structure for the preparation of ab initio QM/MM MD simulations. Each QM/MM model was prepared by deleting the solvent molecules beyond 30 Å from the zinc atom. The resulting QM/MM system had a total of ∼13 000 atoms. The detailed QM/MM partition for all these TLN and HDAC8 models are presented in Figure 1 of the Supporting Information. The QM subsystem, including the zinc and its coordinating ligands, was treated with B3LYP functional using Stuttgart ECP/basis set36 for the zinc atom and 6-31G(d) basis set for all other atoms, which has been previously tested and employed successfully to describe zinc coordination shell.14,21,58-60 The QM/MM boundary was described by improved pseudobond approach.61-64 All other atoms were described by the same molecular mechanical force field used in classical MD simulations. The spherical boundary condition had been applied so that atoms beyond 22 Å from the zinc atom were fixed. The 18 and 12 Å cutoffs were employed for electrostatic and van der Waals interactions, respectively. There was no cutoff for electrostatic interactions between QM and MM regions. The prepared system was first minimized by QM/MM calculations. Finally, more than 20 ps ab initio QM/MM MD simulations were carried out with 1 fs as the time step, and the Beeman algorithm65 was used to integrate the Newton equations of motion. To further check the convergence of our simulations, we also extended the ab initio QM/MM MD simulations to 40 ps for the TLNa system. As shown in Table 1 of the Supporting Information, the difference of average distances and the fluctuation for the zinc coordination shell in TLNa are very similar for different time periods. The Berendsen

340

J. Chem. Theory Comput., Vol. 6, No. 1, 2010

Wu et al.

Figure 2. The representative active site of TLN complex (TLNa) and HDAC8 (model 1). The coordination shell is indicated by a dash line in orange (the most flexible bond is highlighted in purple), and the important hydrogen bonds are shown in green.

thermostat method50 was used to control the system temperature at 300 K. The last 20 ps trajectory, which has achieved the equilibrated temperature (300 K), is adopted as data analysis. To make sure that the fluctuation of the zinc coordination shell is not due to the higher temperature of the QM subsystem, we have monitored it along our QM/ MM MD simulation trajectories. We found that most of the time the temperature of the QM subsystem is not higher than 300 K, and the fluctuation is reasonable. All ab initio QM/ MM calculations were performed with modified Q-Chem66 and Tinker67 programs.

Results and Discussion For TLN, we simulated four enzyme complexes of TLN with different coordinating ligands starting from their respective crystal structures as illustrated in Figure 1a. In the model

Figure 3. The coordination number change during the DFT QM/MM MD simulations.

Thermolysin and HDAC8 Catalytic Zinc Coordination

J. Chem. Theory Comput., Vol. 6, No. 1, 2010 341

Figure 4. The change of the coordination distance to zinc during the DFT QM/MM MD simulations.

TLNa, in which no substrate is bound to the active site, the zinc ion is six-coordinated with a bidentate carboxyl group and two water molecules in the crystal structure. From Table 1, we can see that our MD simulation results not only reproduced the coordination shell very well but also clearly indicated the flexibility and dynamics of the zinc catalytic site. The coordination between Zn2+ and Glu166:OE1 is demonstrated to be most dynamic, and the zinc-carboxylate coordination fluctuates between mono- and bidentate manners as shown in Figures 3 and 4. Interestingly, in the model TLNb, in which the inhibitor is bound to Zn2+ with a thiolate, its zinc coordination is much less flexible which is manifested by the small standard deviation (SD) and is kept to be tetrahedral as in the crystal structure most of the time (see Figure 3). Meanwhile, in models TLNc and TLNd, which have been suggested to mimic the reaction transition state and product, respectively, their zinc active sites are found

to be even more flexible than that in model TLNa. For both models, the zinc-carboxylate coordination is the most flexible and continuously changes between monodentate and bidentate manner dynamically. We can also see that the ligand exchange in zinc coordination is very fast and can occur at the picoseconds scale as illustrated in Figure. 4. For HDAC8, a key distinct feature of its active site is that there are two carboxylate groups coordinated to zinc instead of only one carboxylate in TLN. From structural studies, only monodentate mode of carboxylate coordination has been observed.23-26 Considering the medical and biological significance of HDACs, it would be of particular interest to probe the dynamics of its zinc active site. Here, we have simulated two HDAC8 complexes starting from their respective crystal structures: model 1 is a Y306F mutant in complex with its substrate acetyl-lysine,25 and model 2 is the wildtype HDAC8 binding with its superstar inhibitor SAHA24

342

J. Chem. Theory Comput., Vol. 6, No. 1, 2010

as shown in Figure 1b. From Table 2, we can see that the average distances are consistent with the experimental data, and fluctuations indicate the flexibility of the zinc coordination sphere. Although the zinc active sites in both models are still quite dynamic, the carboxylate coordination is very stable, which is manifested by the small fluctuation, and is always kept in the monodentate mode. As shown in Figures 3 and 4, in the mutant-substrate complex (model 1), its zinc coordination number fluctuates between 4 and 5 and its main flexibility comes from the coordination of the water molecule. In the HDAC8-SAHA complex, its 5-fold zinc coordination is considerably more stable. SAHA is coordinated to the zinc in the bidentate manner most of the time, but it remains to be more flexible than other amino acid ligands, which is indicated by the relatively large standard deviation in Table 2. Thus, our simulations clearly indicate that the dynamics of the zinc active site of HDAC8 is quite different from thermolysin and suggest that HDAC8 is not likely to employ the carboxylate-shift mechanism in its catalytic cycle. As noted in Figure 1, the first coordination shell in TLN is 2His+Glu, and for HDAC8 is 1His+2Asp, so the total charge for the coordination shell is different. This would be an important factor which leads to their different flexibility of zinc-binding sites. Such distinct dynamics in the zinc coordination shell also suggests that the catalytic role of zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of amide bond. In summary, with ab initio QM/MM MD simulations, we have provided direct evidence regarding the inherent flexibility of the catalytic zinc coordination in both TLN and HDAC8 and have observed different coordination modes and fast ligand exchange. For TLN, the coordination of the carboxylate group of Glu166 is found to be most flexible, which can continuously change between monodentate and bidentate manner dynamically. For HDAC8, its coordination to all three amino acid ligands, including two carboxylate groups of Asp, is very stable. Its flexibility mainly comes from a nonamino acid ligand. Such distinct dynamics in their zinc coordination shell suggest that the catalytic role of zinc in TLN and HDAC8 is likely to be different in spite of the fact that both catalyze the hydrolysis of amide bond. Furthermore, our work demonstrates the feasibility and applicability of Born-Oppenheimer ab initio QM/MM MD simulations in simulating metalloenzymes and sets the stage for more detailed understanding of catalysis and inhibition in TLN and HDAC8. Acknowledgment. This work was supported by NIH (R01-GM079223), NSF (CHE-CAREER-0448156), and the China Scholarship Council. We thank NYU-ITS and NCSA for providing computational resources. Supporting Information Available: The illustration of QM/MM partition for selected model systems and the data for 40 ps ab initio QM/MM MD results. This information is available free of charge via the Internet at http://pubs.acs.org

Wu et al. (2) Odell, B. L. Nutr. ReV. 1992, 50, 48–50. (3) Lipscomb, W. N.; Strater, N. Chem. ReV. 1996, 96, 2375– 2433. (4) Christianson, D. W.; Cox, J. D. Annu. ReV. Biochem. 1999, 68, 33–57. (5) Parkin, G. Chem. ReV. 2004, 104, 699–767. (6) Anzellotti, A. I.; Farrell, N. P. Chem. Soc. ReV. 2008, 37, 1629–1651. (7) Karlin, S.; Zhu, Z. Y. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 14231–14236. (8) Lee, Y. M.; Lim, C. J. Mol. Biol. 2008, 379, 545–553. (9) McCall, K. A.; Huang, C. C.; Fierke, C. A. J. Nutr. 2000, 130, 1437S–1446S. (10) Sigel, H.; Martin, R. B. Chem. Soc. ReV. 1994, 23, 83–91. (11) Parkin, G. Chem. Commun. 2000, 20, 1971–1985. (12) Tamames, B.; Sousa, S. F.; Tamames, J.; Fernandes, P. A.; Ramos, M. J. Proteins: Struct., Funct., Bioinf. 2007, 69, 466–475. (13) Robert, V.; Lemercier, G. J. Am. Chem. Soc. 2006, 128, 1183–1187. (14) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. J. Am. Chem. Soc. 2007, 129, 1378–1385. (15) Matthews, B. W. Acc. Chem. Res. 1988, 21, 333–340. (16) Holland, D. R.; Hausrath, A. C.; Juers, D.; Matthews, B. W. Protein Sci. 1995, 4, 1955–1965. (17) Roderick, S. L.; Fourniezaluski, M. C.; Roques, B. P.; Matthews, B. W. Biochemistry 1989, 28, 1493–1497. (18) Tronrud, D. E.; Monzingo, A. F.; Matthews, B. W. Eur. J. Biochem. 1986, 157, 261–268. (19) Holden, H. M.; Matthews, B. W. J. Biol. Chem. 1988, 263, 3256–3260. (20) Kimura, E. Acc. Chem. Res. 2001, 34, 171–179. (21) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. Biophys. J. 2005, 88, 483–494. (22) Szeto, M. W. Y.; Mujika, J. I.; Zurek, J.; Mulholland, A. J.; Harvey, J. N. J. Mol. Struct.: THEOCHEM 2009, 898, 106– 114. (23) Vannini, A.; Volpari, C.; Filocamo, G.; Casavola, E. C.; Brunetti, M.; Renzoni, D.; Chakravarty, P.; Paolini, C.; De Francesco, R.; Gallinari, P.; Steinkuhler, C.; Di Marco, S. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 15064–15069. (24) Somoza, J. R.; Skene, R. J.; Katz, B. A.; Mol, C.; Ho, J. D.; Jennings, A. J.; Luong, C.; Arvai, A.; Buggy, J. J.; Chi, E.; Tang, J.; Sang, B. C.; Verner, E.; Wynands, R.; Leahy, E. M.; Dougan, D. R.; Snell, G.; Navre, M.; Knuth, M. W.; Swanson, R. V.; McRee, D. E.; Tari, L. W. Structure 2004, 12, 1325– 1334. (25) Vannini, A.; Volpari, C.; Gallinari, P.; Jones, P.; Mattu, M.; Carfi, A.; De Francesco, R.; Steinkuhler, C.; Di Marco, S. EMBO Rep. 2007, 8, 879–884. (26) Dowling, D. P.; Gantt, S. L.; Gattis, S. G.; Fierke, C. A.; Christianson, D. W. Biochemistry 2008, 47, 13554–13563.

References

(27) Mann, B. S.; Johnson, J. R.; Cohen, M. H.; Justice, R.; Pazdur, R. Oncologist 2007, 12, 1247–1252.

(1) Andreini, C.; Banci, L.; Bertini, I.; Rosato, A. J. Proteome Res. 2006, 5, 196–201.

(28) Marks, P. A.; Breslow, R. Nat. Biotechnol. 2007, 25, 84– 90.

Thermolysin and HDAC8 Catalytic Zinc Coordination

J. Chem. Theory Comput., Vol. 6, No. 1, 2010 343

(29) Park, H.; Brothers, E. N.; Merz, K. M. J. Am. Chem. Soc. 2005, 127, 4232–4241.

P. A. AMBER 10 University of California: San Francisco, CA, 2008.

(30) Wong, K. Y.; Gao, J. Biochemistry 2007, 46, 13352–13369.

(52) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179– 5197.

(31) Elstner, M.; Cui, Q.; Munih, P.; Kaxiras, E.; Frauenheim, T.; Karplus, M. J. Comput. Chem. 2003, 24, 565–581. (32) Lopez-Canut, V.; Marti, S.; Bertran, J.; Moliner, V.; Tunon, I. J. Phys. Chem. B 2009, 113, 7816–7824. (33) Brauer, M.; Kunert, M.; Dinjus, E.; Klussmann, M.; Doring, M.; Gorls, H.; Anders, E. J. Mol. Struct.: THEOCHEM 2000, 505, 289–301. (34) Lopez, X.; York, D. M. Theor. Chem. Acc. 2003, 109, 149– 159. (35) Xu, D. G.; Guo, H. J. Am. Chem. Soc. 2009, 131, 9780– 9788. (36) Dolg, M.; Wedig, U.; Stoll, H.; Preuss, H. J. Chem. Phys. 1987, 86, 866–872. (37) Hu, P.; Wang, S.; Zhang, Y. J. Am. Chem. Soc. 2008, 130, 3806–3813. (38) Hu, P.; Wang, S.; Zhang, Y. J. Am. Chem. Soc. 2008, 130, 16721–16728. (39) Wang, S. L.; Hu, P.; Zhang, Y. K. J. Phys. Chem. B 2007, 111, 3758–3764.

(53) Wang, J. M.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049–1074. (54) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins: Struct., Funct., Bioinf. 2006, 65, 712–725. (55) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (56) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157–1174. (57) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (58) Corminboeuf, C.; Hu, P.; Tuckerman, M. E.; Zhang, Y. K. J. Am. Chem. Soc. 2006, 128, 4530–4531. (59) Xiao, C. Y.; Zhang, Y. K. J. Phys. Chem. B 2007, 111, 6229– 6235. (60) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. Proteins: Struct., Funct., Bioinf. 2007, 66, 205–218.

(40) Ke, Z. H.; Zhou, Y. Z.; Hu, P.; Wang, S. L.; Xie, D. Q.; Zhang, Y. K. J. Phys. Chem. B 2009, 113, 12750–12758.

(61) Zhang, Y. K. Theor. Chem. Acc. 2006, 116, 43–50.

(41) Guex, N.; Peitsch, M. C. Electrophoresis 1997, 18, 2714– 2723.

(62) Zhang, Y. K.; Lee, T. S.; Yang, W. T. J. Chem. Phys. 1999, 110, 46–54.

(42) Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. Nucleic Acids Res. 2005, 33, W368–W371.

(63) Zhang, Y. K.; Liu, H. Y.; Yang, W. T. J. Chem. Phys. 2000, 112, 3483–3492.

(43) Antonczak, S.; Monard, G.; Ruiz-Lopez, M. F.; Rivail, J. L. J. Am. Chem. Soc. 1998, 120, 8825–8833.

(64) Zhang, Y. K. J. Chem. Phys. 2005, 122, 024114.

(44) Blumberger, J.; Lamoureux, G.; Klein, M. L. J. Chem. Theory Comput. 2007, 3, 1837–1850. (45) Stote, R. H.; Karplus, M. Proteins: Struct., Funct., Genet. 1995, 23, 12–31. (46) Liang, J. Y.; Lipscomb, W. N. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 3675–3679. (47) Pang, Y. P. Proteins: Struct., Funct., Genet. 2001, 45, 183– 189. (48) Yan, C. L.; Xiu, Z. L.; Li, X. H.; Li, S. M.; Hao, C.; Teng, H. Proteins: Struct., Funct., Bioinf. 2008, 73, 134–149. (49) Karjiban, R. A.; Rahman, M. B. A.; Basri, M.; Salleh, A. B.; Jacobs, D.; Wahab, H. A. Protein J. 2009, 28, 14–23. (50) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (51) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W.; Merz, K. M.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossva´ry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V. Kollman,

(65) Beeman, D. J. Comput. Phys. 1976, 20, 130–139. (66) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.; DiStasio, R. A.; Lochan, R. C.; Wang, T.; Beran, G. J. O.; Besley, N. A.; Herbert, J. M.; Lin, C. Y.; Van Voorhis, T.; Chien, S. H.; Sodt, A.; Steele, R. P.; Rassolov, V. A.; Maslen, P. E.; Korambath, P. P.; Adamson, R. D.; Austin, B.; Baker, J.; Byrd, E. F. C.; Dachsel, H.; Doerksen, R. J.; Dreuw, A.; Dunietz, B. D.; Dutoi, A. D.; Furlani, T. R.; Gwaltney, S. R.; Heyden, A.; Hirata, S.; Hsu, C. P.; Kedziora, G.; Khalliulin, R. Z.; Klunzinger, P.; Lee, A. M.; Lee, M. S.; Liang, W.; Lotan, I.; Nair, N.; Peters, B.; Proynov, E. I.; Pieniazek, P. A.; Rhee, Y. M.; Ritchie, J.; Rosta, E.; Sherrill, C. D.; Simmonett, A. C.; Subotnik, J. E.; Woodcock, H. L.; Zhang, W.; Bell, A. T.; Chakraborty, A. K.; Chipman, D. M.; Keil, F. J.; Warshel, A.; Hehre, W. J.; Schaefer, H. F.; Kong, J.; Krylov, A. I.; Gill, P. M. W.; HeadGordon, M. Q-Chem, version 3.0; Q-Chem, Inc.: Pittsburgh, PA, 2006. (67) Ponder, J. W. TINKER, Software Tools for Molecular Design, version 4.2; 2004. (68) Maret, W.; Li, Y. Chem. ReV. 2009, 109, 4682–4707. CT9005322