J. Phys. Chem. B 2007, 111, 12909-12915
12909
Comparison of Different Quantum Mechanical/Molecular Mechanics Boundary Treatments in the Reaction of the Hepatitis C Virus NS3 Protease with the NS5A/5B Substrate Alejandro Rodrı´guez, Carolina Oliva,* and Miguel Gonza´ lez* Departament de Quı´mica Fı´sica i Centre de Recerca en Quı´mica Teo` rica, UniVersitat de Barcelona i Parc Cientı´fic de Barcelona, C/Martı´ i Franque` s 1, 08028 Barcelona, Spain
Marc van der Kamp and Adrian J. Mulholland* Centre for Computational Chemistry, School of Chemistry, UniVersity of Bristol, Bristol BS8 1TS, U.K. ReceiVed: June 5, 2007; In Final Form: August 6, 2007
The link atom (LA) and the generalized hybrid orbital (GHO) quantum mechanical/molecular mechanics (QM/MM) boundary treatment methods are compared, in the context of the acylation process (the ratelimiting step) involving the NS3/NS4A HCV serine protease and its NS5A/5B natural substrate. The potential energy surface was calculated, and the free energy along the selected reaction coordinate was obtained from umbrella sampling molecular dynamics simulations, at the AM1/CHARMM27 level. The LA and GHO methods, when properly applied, lead to similar chemical behavior, although the agreement is not quantitative. The choice of QM/MM partitioning is dictated to some extent by the nature of the two different methods, and this influences the results. The free energy profiles obtained by umbrella sampling suggest that the GHO approach is better suited for this system, because it provides a consistent description of the reaction in both the forward and backward directions. This is probably a consequence of the different QM/MM partitioning required by the two different methods (i.e., different numbers of atoms have to be included in the QM region). This finding is therefore likely to be system dependent, so careful testing should be considered for each enzyme application.
Introduction Knowledge of the detailed nature of the chemical reactivity in large systems, e.g., reactions that take place in enzymes or in solution, has become one of the main issues in theoretical chemistry. The size makes it computationally unfeasible to study this kind of systems with conventional quantum mechanical (QM) methods, except by considering reduced models. The changes in the electronic structure and bonding involved in a chemical reaction also preclude the use of standard molecular mechanics (MM) methods. Thus, the development of alternative treatments to deal with these large, and often complex, systems is an active field, despite the increasing computational capability available. Among the most popular and powerful approaches are those based on hybrid methods, where the system is divided into two subsystems, namely, the subsystem directly involved in reactivity (treated with a QM method) and the rest of the system (treated with an MM method).1-3 These QM/MM methods can be classified in several ways, e.g., depending on the general hybrid energy expression4 or on the kind of coupling between both subsystems.5 It is now possible to perform QM/MM calculations (when using high levels of electronic structure description, e.g., the Coupled Cluster/MM) that can give chemical accuracy for activation barriers of enzyme-catalyzed reactions.6 * Corresponding authors. Tel.: +44 1223 228664 (C.O.); +34 93 402 12 37 (M.G.); +44 (0)117 928 9097 (A.J.M.). Fax: +44 1223 228501 (C.O.); +34 93 402 12 31 (M.G.); +44 (0)117 925 1295 (A.J.M.). E-mail:
[email protected] (C.O.);
[email protected] (M.G.);
[email protected] (A.J.M.).
Whatever the hybrid method used, a nontrivial problem appears when the division between the QM and MM parts cuts a covalent bond. Then, the standard QM/MM interaction models are not enough to describe the strong electronic structure change produced. In fact, if the charge of the QM part is conserved, this implies a change from a closed shell structure to an open shell one and, due to this, a special treatment of the boundary has to be applied (see, e.g., ref 7). The different approaches to the boundary treatment problem can be grouped in two classes. In the first, an atom (link atom, connection atom, or a pseudoatom) with a free valence available is added to saturate the cut bond.8-11 In the second class, the treatments try to reproduce directly the characteristics of the cut bond by constraining the SCF (self-consistent field) solution,12-15 with the purpose of avoiding the introduction of extra geometric degrees of freedom implied in the first class of methods. There is no consensus about which is the best way to deal with this problem, it being generally accepted that all approaches have advantages and drawbacks. Several papers studied the effect of the treatment in the boundary,16-18 but mainly comparing with ab initio results in small systems. As pointed out by Lin and Truhlar recently,3 it is very hard to test schemes designed for processes that occur in condensed phase by only carrying out calculations on small molecular models in the gas phase. Thus, condensed phase test cases are needed in order to understand the complexities inherent when working with each scheme.17 Here, we have focused on two widely used methods for the QM/MM boundary treatment, each of them being representative
10.1021/jp0743469 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/13/2007
12910 J. Phys. Chem. B, Vol. 111, No. 44, 2007
Rodrı´guez et al.
Figure 1. MM charge modification in one example employing the HQ link atom type. See text.
of one of the two classes mentioned above (since the two classes include quite different methods, the selection of a representative one is necessarily criterion dependent). We have compared the behavior of the link atom (LA) and generalized hybrid orbital (GHO) approaches for QM/MM modeling of a real enzyme reaction: we believe that it is important to perform such tests on enzyme systems, rather than solely on small model systems. In the GHO method, a set of hybrid orbitals is placed on the boundary atom (the QM atom covalently bonded to the MM region, always an sp3 carbon atom) between the QM and MM fragments, and one of the hybrid orbitals participates in the SCF calculation made for the atoms in the QM region.14 There are several types of LA schemes.8,17 All of them have in common the use of a new QM atom (that does not exist in the real system, usually a hydrogen atom), placed along the cut QM/MM bond and which adds one electron to saturate the free valence. In the LA implementation (“QQ” type)9 of Field et al., the electrostatic QM/MM interaction terms were set to zero when the QM atom is the link atom. However, it was stated that the different treatment of subsets inside the QM region led to the creation of an artificial charge distribution in this region, and this, while modifying the polarities, can adversely affect the computed energies.10 Thus, here we used the “HQ” link atom type,10 where the link atoms are allowed to interact with the MM part and the charges of the surrounding MM atoms are modified in a proper way for the QM/MM calculation (see an example in Figure 1 to clarify this point). The charge of the MM neighbor atom to the link atom is set to zero (R carbon in the example shown in Figure 1), and the rest of the electrostatic group MM charges are modified conveniently to maintain the charge of the group in the CHARMM22 force field.19 In this force field, the electrostatic group, where the R carbon is located, is composed by the R carbon, R hydrogen bonded to the R carbon, the backbone nitrogen atom, and the hydrogen atom bonded to this nitrogen atom. Due to the need for condensed phase test cases, here we have investigated at the QM/MM level the first step (acylation process; Figure 2) of the hydrolysis of the NS5A/5B junction catalyzed by the hepatitis C virus NS3/NS4A serine protease, which corresponds to the rate-limiting step.20 The effect of the QM/MM boundary treatment on the potential energy surface (PES) and the relative free energy of the system along a selected reaction coordinate has been explored.
Figure 2. Acylation process (rate-limiting step) for the HCV NS3/ NS4A protease catalytic reaction mechanism.
Computational Methods The initial geometry considered was a 2.4 Å resolution structure of a NS3/NS4A protease-inhibitor complex (PDB code 1DXP21), while to model the NS5A/5B substrate,22 we overlapped the catalytic triad (His-57, Asp-81, and Ser-139) and Arg-155 in the above-mentioned structure for NS3, with the corresponding residues (His-57, Asp-102, Ser-195, and Ser-214) in the crystallographic structure of chymotrypsin with Eglin C at 2.0 Å resolution (PDB code 1ACB23), and the residues of the ligand were mutated to the corresponding ones in the NS5A/ 5B substrate. A “cationic dummy atom” model24 was used to model the Zn2+ present in the reactant structure and far from the active site, while the parameters of the AMBER25 force field were employed for the three deprotonated cysteine residues bonded to the zinc cation. The remaining residues were modeled with the CHARMM2219 all-atom parameter set. Finally, a sphere of pre-equilibrated CHARMM TIP3P26 water molecules of radius 30 Å, centered on the center of mass of the whole system, was used to solvate the enzyme-substrate complex. All computations have been done with the c27b2 and c29b1 versions of the CHARMM program.27 A more detailed description of the setup of the system can be found in ref 20. This structure was divided into two subsystems (QM and MM regions). The QM region includes part of the substrate and the side chains of the enzyme residues, and three different partition schemes were applied (Figure 3). In all calculations, the QM region was treated using the semiempirical method AM1, and the wave function is polarized by the MM environment (electrostatic embedding). The criteria for the selection of the QM region were the following: (a) in the GHO case, the division was made ensuring that all GHO atoms correspond to sp3 carbon atoms, since it is a method requirement; (b) in the first LA partition scheme, the objective was to obtain the maximum similarity with the corresponding GHO partition, so the LA partition was done along the same bonds modeled by the active hybrid orbitals of the GHO partition; (c) the second LA partition scheme allowed us to test of the influence of
Comparison of Different QM/MM Boundary Treatments
Figure 3. Different partition schemes used in this study. (A) GHO, (B) first LA, and (C) second LA partition schemes. Atom labeling is shown in (A).
changing the location of the link atoms (the GHO approach is not so flexible, because the GHO atom must currently be an sp3 carbon). In particular, the link atom located along the bond between the backbone nitrogen and the R carbon of the Met residue situated in the site P2′ is a candidate to generate problems, due to the difference of polarities between these two atoms. Because of this, two new link atoms were added to ensure that all the cut bonds were between carbons. Thus, 53 atoms are included in the QM region in the GHO and first LA partition schemes, while in the second LA partition scheme the number of QM atoms was 56. In all cases a total QM charge of -1 is used. The reactant structures were obtained by energy minimizing the initial structure with the steepest descent (SD) and the adopted basis Newton-Raphson27,28 (ABNR) methods, until the root-mean-square value of the components of the gradient vector in a Cartesian coordinate system (GRMS) fell below 0.001 kcal mol-1 Å-1. In all cases, the 30 Å pseudosphere of CHARMM
J. Phys. Chem. B, Vol. 111, No. 44, 2007 12911 TIP3P water molecules was allowed to relax in the minimization under a deformable spherical boundary potential29 to prevent water from escaping into the vacuum. The next step was to study the PES obtained from the three methods. To this end, two distinguished reaction coordinates (DCs) were defined to describe the acylation process, as in ref 20: (a) the proton transfer between Ser-139 and His-57 (difference between the Os-Hs and Hs-N distances; see Figure 3 for atom naming); (b) the nucleophilic attack of Ser139 on the carbonyl carbon of the substrate (Os-Cp distance). The PES calculations were carried out by restraining independently the DCs to 0.05 Å separated values from reactants to products with a harmonic potential (force constant of 2000 kcal mol-1 Å-2). The energy was minimized for each selected value of the DC by means of the ABNR method until the GRMS fell below 1 × 10-2 kcal mol-1 Å-1. The results of these calculations are shown in Figure 4 and will be considered in detail in the Results and Discussion section. These methods were also employed to obtain the intermediate (I1) and barrier (B1) structures, the first one being energy minimized from the nearest structure until GRMS e 10-3 kcal mol-1 Å-1; and the barrier state (energy maximum along the minimum energy path connecting reactants and the intermediate using the selected coordinates) was further refined by means of the TRAVEL30 module in CHARMM, until the GRMS fell below 0.01 kcal mol-1 Å-1. After analysis of the PES for the acylation process (ratelimiting step of the reaction), the influence of the temperature on this step was analyzed, and the activation free energy was calculated at 300 K, considering both the GHO and the LA (with the second partition scheme) approaches. To this end, an umbrella sampling31 study using SBMD32,33 (stochastic boundary molecular dynamics) was performed, with a reaction region (treated with Newtonian dynamics) of 26 Å radius and a buffer region (treated with Langevin dynamics) from 26 to 30 Å. From the reactants and the intermediate structure, the system was heated for 1.0 ps at 300 K, equilibrated at this temperature for 2.0 ps, and allowed to evolve freely for 10.0 ps. The definition of the reaction coordinate as a combination of distances has been successfully applied for umbrella sampling simulations,34,35 and in this case it has been defined as ζ ) d(OS-CP) + d(NHS) - d(OS-HS), since the variation along the coordinates d(OS-CP) and d(OS-HS) - d(N-HS) described the acylation process. The comparison of the barrier heights obtained scanning this coordinate with those obtained from the PES (35.8 vs 35.6 kcal mol-1, respectively, in the GHO case and 42.2 kcal mol-1 in both cases in the LA one) suggests the goodness of this distinguished coordinate. This coordinate was sampled with the umbrella sampling technique in molecular dynamics simulations, in the interval between 0.0 and 5.0 Å, by restraining it using a harmonic potential with a force constant of 100 kcal mol-1 Å-2 to ensure the surmounting of the barrier. The step of the sampling (window range) was equal to 0.1 Å. For every window in the reaction coordinate, 10 ps of SBMD was carried out, and the WHAM36 (weighted histogram analysis method) technique was used to collate the data from all simulations and calculate the relative free energy. Results and Discussion The analysis of the PES showed a similar general behavior for the three partition schemes, and acylation occurs in a concerted way (Figure 4). The results are in agreement with those obtained for other proteases.37-40 The nonrealistic height of the barriers comes from the trend of AM1 to overestimate
12912 J. Phys. Chem. B, Vol. 111, No. 44, 2007
Rodrı´guez et al. TABLE 1: Summary of the Main Free Energy Differences Results (∆G0) and Selected Energies (∆E) on the Potential Energy Surface (in kcal mol-1 Taken from Reactants)a GHO forward ξ
∆G
second LA partition
backward ξ
∆G
PES ξ
∆E
forward ξ
∆G
backward ξ
∆G
PES ξ
∆E
B 1.2 42.3 1.3 38.2 1.56 35.6 1.5 49.7 1.5 46.1 1.53 42.2 I 0.45 35.0 0.5 30.1 0.55 24.3 0.95 45.2 0.4 34.6 0.52 31.14 a “Forward” and “backward” refer to the ∆G calculation, and “PES” refers to the energy of the stationary points on the PES. The ∆G and ∆E values are taken with respect to reactants.
Figure 4. Equipotential contour plots of the two-dimensional representation of the PES for the acylation process computed with the three partition schemes. (A) GHO, (B) first LA, and (C) second LA.
them. In ref 20 a correction of the energies using MP2/6-31G* calculations on the AM1/MM obtained geometries gives more reliable values (about 10 kcal mol-1). In both cases, the refined barrier structures are really similar to those obtained by restraining the reaction coordinate. However, there are some differences that should be mentioned. When comparing both LA partition schemes with GHO, there is a barrier displacement, and the barrier height increases in about 6 kcal mol-1 (41.6
kcal mol-1 in the first LA case and 42.2 kcal mol-1 in the second one vs 35.6 kcal mol-1 in the GHO case). Moreover, the first PES of the LA partition scheme presents a significant roughness near the intermediate (a new minimum appears to be formed). An explanation of the last point could come from a detailed structural analysis of the intermediates found on the PES (see Figure 5). From these figures, it is evident that the umbrella motion in the Np nitrogen has not occurred yet in the GHO treatment, while in the first LA case it has already occurred (for this reason, the intermediate in Figure 4B has been denoted as I and not as I1). From our previous study on this system,20 it is known that in the GHO case this motion appears to happen after finishing the acylation process and before the peptide bond breaking, which is the next step of the catalytic mechanism. However, it is widely known that the AM1 method is not always reliable for modeling nitrogen inversion, so this difference could be an artifact of the method. Analysis of the structure of the intermediate from the PES with the second LA partition scheme confirmed that in this case the umbrella motion in the Np nitrogen has not occurred. This and the fact that the PES wrinkling around the intermediate has disappeared (a new roughness has appeared, but it is far from the stationary points), suggests that it could be an artifact coming from a bad choice of the QM/MM partition. Thus, from now on, we will discuss only with the second LA partition. In such a complex system, it is difficult to distinguish small changes in PES for reaction (e.g., due to QM/MM partitioning schemes) from small changes in conformation (e.g., preferences for different local minima). Comparisons based on single PESs (as opposed to large numbers of conformations) should be made with caution. Nevertheless, it appeared that this first choice of QM/MM partition gave less satisfactory results. The results of the WHAM analysis of the umbrella sampling calculations starting from reactants (forward reaction) and from the intermediate (backward reaction) are shown in Figure 6 and summarized in Table 1 for GHO and the second LA partition. The locations for the barrier and the intermediate ensembles have been denoted as B and I, respectively. The free energy values for the barrier are bigger than the potential energy ones, as expected. As discussed in ref 20, the process corresponds to an intermolecular bond formation, so the entropic term is expected to be negative, and consequently, the free energy will increase with temperature. From an exploration of the trajectory geometries, it can be concluded that the difference found between the behavior in the intermediate region of the free energy profiles of the reaction for the forward and backward pathways for the second LA partition (the forward pathway does not describe the intermediate structure which is the starting point for the backward one) arises from an internal rotation in the Asp-81 residue of the catalytic triad in the reactant structure (Figure 7). Due to this internal
Comparison of Different QM/MM Boundary Treatments
J. Phys. Chem. B, Vol. 111, No. 44, 2007 12913
Figure 5. Comparison of the intermediate structures found for the GHO (A) and the first LA (B) partition schemes.
Figure 7. Illustration of the internal rotation of the Asp-81 residue in the reactant structure.
Figure 6. Free energy profiles from the WHAM analysis of the umbrella sampling calculations for the GHO (A) and the second LA (B) partition schemes at the AM1/CHARMM level.
rotation, the hydrogen bond between the Asp-81 and the His57 is broken, and the stabilization of the charged His-57 in the intermediate (and in the barrier state) does not occur. Since the tetrahedral intermediate has been described as necessary for the evolution of the catalytic mechanism in other serine proteases,41 the destabilization of this intermediate, when using the second LA partition, seems to be a drawback for the study of this
particular system with this partition method. The internal rotation in the Asp-81 residue is considered in more detail below. To study the origin of this difference, energy scans along the dihedral angle that describes the rotation around the CRCβ bond of this Asp-81 residue (Figure 7) were carried out (Figure 8A). It was restrained with a harmonic potential (force constant of 2000 kcal mol-1 Å-2) to a series of values separated by 5° and energy minimized by means of the ABNR method (GRMS e 10-2 kcal mol-1 Å-1). From Figure 8A, it can be seen that there is an energetic basis for the different behavior of the system as a function of the dihedral angle, and this depends on the QM/MM partitioning method employed. The second LA partition method gives a smaller rotational barrier when compared with the GHO result. This finding is consistent with the fact that this internal rotation becomes easier in the former case and, thus, explains the results of the umbrella sampling dynamics. It is important to make clear that the QM systems are not the same for the two models, because of the need for different partitioning for the GHO and LA methods. The origin of the observed difference may well be in the
12914 J. Phys. Chem. B, Vol. 111, No. 44, 2007
Rodrı´guez et al.
Figure 8. Energy profiles for the internal rotation of the aspartate residue: (A) scan for Asp-81 in the reactant structure with the GHO and the second LA partition schemes at the AM1/CHARMM level; (B) scan in the model structure shown with pure MM (CHARMM27), pure AM1, and hybrid GHO and second LA partition schemes.
different QM (and MM) systems, rather than because of the partitioning schemes themselves. To achieve a deeper understanding of the behavior indicated above, regarding the internal rotation of Asp-81, the model system showed as inserto in Figure 8B was also studied (an aspartate residue with an acetylated N-terminus and a methylated C-terminus). The internal rotation was investigated in the model employing the MM, AM1, and the hybrid AM1/MM methods using both the GHO approach (over the CR atom) and the LA partition approach (we have considered a single LA option, the second model, where the HQ atom is situated along the CRCβ bond) to saturate the QM part, the parameters for the angle scanning being identical to the previous ones. The results are shown in Figure 8B, and representative structures for the minima found are given in Figure 9. The second LA approach behaves quite similarly to the pure MM method (CHARMM), whereas the GHO energy curve is more similar to that for the pure AM1 method. These findings and, in particular, the backward behavior found when studying the model with respect to what was obtained for the full system, confirm the crucial importance of the division scheme, even in a system as simple as the aspartate
Figure 9. Minima found (5°, 190°, and 310°) for the internal rotation of the aspartate residue model using the MM (CHARMM27) Hamiltonian.
model analyzed. This is, of course, not trivial for complex systems. The results here show that different partitioning schemes may enforce different choices of QM region, and this choice may have a significant effect on the conformational or energetic behavior of the system. Different types of QM/MM partitioning method may require different sizes of QM region, and different locations for the QM/MM boundary. This may lead to significant differences in behavior, not caused by the QM/MM partition itself. It is worth pointing out that a larger QM region is not always desirable; particularly when lowintermediate levels of QM theory are used, structure (e.g.,
Comparison of Different QM/MM Boundary Treatments protein) and interactions may in fact be better treated by a good MM force field. Summary and Conclusions A comparison of the behavior of the LA and GHO QM/MM boundary methods for the acylation process (rate-limiting step) of the NS3/NS4A HCV serine protease catalytic mechanism with one of its main natural substrates (NS5A/5B) has been performed, with the purpose to test these methods on a real application. To this end, calculations on the PES and dynamic simulations with different partition schemes have been carried out, considering also that different sizes of the QM region could be adapted better to each of the partition schemes considered. A test has also been carried out on a model system, to interpret the results derived from umbrella sampling calculations. We conclude that the results for this reaction are dependent on the QM/MM partitioning method used, so this aspect should be carefully tested in QM/MM investigations. Moreover, the inclusion of temperature, with the increased conformational sampling that it implies, can increase the observed differences. The LA method is more flexible than the GHO method (since the GHO atom must be an sp3 carbon), but it also introduces some uncertainty in the link atom location. To maintain the polarity of the cut bonds seems to be crucial, so it is recommended to cut only C-C bonds, since the electronegativities of C and H are similar enough generally to model a C-C bond with a C-H one. However, both the LA and GHO partition methods, when properly applied, lead to similar behavior of this system (i.e., nonsequential mechanism for the acylation process, increase of the free energy barrier with temperature, etc.). Hence, although it is hard to conclude which is best suited for this enzyme system, the umbrella sampling results suggest that the GHO approach is better in this case. For other systems, however, things may be different. Acknowledgment. This work was supported by the Spanish Ministry of Education and Science (Project CTQ2005-09334C02-01). Thanks are also given to the “Generalitat” (Autonomous Government) of Catalonia (Project 2005SGR 00175), and A.R. thanks the Spanish Ministry of Science and Technology for a predoctoral research Grant (this work is a partial fulfillment of his Ph.D. Thesis). Finally, the authors are grateful to the “Centre de Supercomputacio´ de Catalunya” (CESCA) for providing part of the computer time. A.J.M. thanks the IBM High Performance Computing Life Sciences Outreach Program, BBSRC, and EPSRC. M.v.d.K. is funded by the BRISENZ Marie Curie Early Stage Training Centre: he and A.J.M. thank the European Commission for this support. References and Notes (1) Mulholland, A. J. Drug DiscoVery Today 2005, 10, 1393. (2) Friesner, R. A.; Guallar, V. Annu. ReV. Phys. Chem. 2005, 56, 389. (3) Lin, H.; Truhlar, D. G. Theor. Chem. Acc. 2007, 117, 185. (4) Sherwood, P. Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC Series; Forschungszentrum Juelich, Germany, 2000; Vol. 3, p 285.
J. Phys. Chem. B, Vol. 111, No. 44, 2007 12915 (5) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (6) Claeyssens, F.; Harvey, J. N.; Manby, F. R.; Mata, R. A.; Mulholland, A. J.; Ranaghan, K. E.; Schu¨tz, M.; Thiel, S.; Thiel, W.; Werner, H.-J. Angew. Chem., Int. Ed. 2006, 45, 6856. (7) Wang, W.; Donini, O.; Reyes, C. M.; Kollman, P. A. Annu. ReV. Biophys. Biomol. Struct. 2001, 30, 211. (8) Antes, I.; Thiel, W. ACS Symp. Ser. 1998, 712, 50. (9) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (10) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M. J. Phys. Chem. A 2000, 104, 1720. (11) DiLabio, G. A.; Hurley, M. M.; Christiansen, P. A. J. Chem. Phys. 2002, 116, 9578. (12) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (13) The´ry, V.; Rinaldi, D.; Rivail, J. L.; Maigret, B.; Frenczy, B. J. Comput. Chem. 1994, 15, 269. (14) Gao, J.; Amara, P.; Alhambra, C.; Field, M. J. J. Phys. Chem. A 1998, 12, 4714. (15) Murphy, R.; Philipp, D.; Friesner, R. Chem. Phys. Lett. 2000, 321, 113. (16) Lennartz, C.; Scha1fer, A.; Terstegen, F.; Thiel, W. J. Phys. Chem. B 2002, 106, 1758. (17) Konig, P. H.; Hoffmann, M.; Frauenheim, T.; Cui, Q. J. Phys. Chem. B 2005, 109, 9082. (18) Nemukhin, A. V.; Grigorenko, B. L.; Topol, I. A.; Burt, S. K. J. Comput. Chem. 2003, 24, 1410. (19) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (20) Oliva, C.; Rodrı´guez, A.; Gonza´lez, M.; Yang, W. Proteins: Struct., Funct., Bioinf. 2007, 66, 444. (21) Di Marco, S.; Rizzi, M.; Volpari, C.; Walsh, M.; Narjes, F.; Colarusso, S.; De Francesco, R.; Matassa, V. G.; Sollazzo, M. J. Biol. Chem. 2000, 275, 7152. (22) Kwong, A. D.; Kim, J. L.; Rao, G.; Lipovsek, D.; Raybuck, S. A. AntiViral Res. 1998, 40, 1. (23) Frigerio, F.; Coda, A.; Pugliese, L.; Lionetti, C.; Menegatti, E.; Amiconi, G.; Schnebli, H. P.; Ascenzi, P.; Bolognesi, M. J. Mol. Biol. 1992, 225, 107. (24) Pang, Y.; Xu, K.; El Yazal, J.; Prendergast, F. G. Protein Sci. 2000, 9, 1857. (25) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E., III; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Comput. Phys. Commun. 1995, 91, 1. (26) Chambers, C. C.; Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 16385. (27) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; Status, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (28) Chu, J. H.; Trout, B. L.; Brooks, B. R. J. Chem. Phys. 2003, 119, 12708. (29) Brooks, C. L.; Karplus, M. J. Chem. Phys. 1983, 79, 6312. (30) Fischer, S.; Karplus, M. Chem. Phys. Lett. 1992, 194, 252. (31) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (32) Brooks, C. L.; Karplus, M. J. Mol. Biol. 1989, 208, 159. (33) Ridder, L.; Rietjens, I. M. C. M.; Vervoort, J.; Mulholland, A. J. J. Am. Chem. Soc. 2002, 124, 9926. (34) Devi-Kesavan, L. S.; Gao, J. J. Am. Chem. Soc. 2003, 125, 1532. (35) Alhambra, C.; Wu, L.; Zhang, Z. Y.; Gao, J. J. Am. Chem. Soc. 1998, 120, 3858. (36) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1992, 13, 1011. (37) Hedstrom, L. Chem. ReV. 2002, 102, 4501. (38) Warshel, A.; Naray-Szabo, G.; Sussman, F.; Hwang, J.-K. Biochemistry 1999, 28, 3629. (39) Ishida, T.; Kato, S. J. Am. Chem. Soc. 2004, 126, 7111. (40) Ishida, T.; Kato, S. J. Am. Chem. Soc. 2003, 125, 12035. (41) Radisky, E. S.; Lee, J. M.; Lu, C. K.; Koshland, D. E., Jr. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 6835.
