2212
J. Phys. Chem. B 2010, 114, 2212–2218
VB/MM Protein Landscapes: A Study of the SN2 Reaction in Haloalkane Dehalogenase† Avital Sharir-Ivry, Tamar Shnerb, Marek Sˇtrajbl, and Avital Shurki*,‡ Department of Medicinal Chemistry and Natural Product, The Institute of Drug Research, The Lise Meitner-MinerVa Center for Computational Quantum Chemistry, The Hebrew UniVersity of Jerusalem, 91120 Jerusalem, Israel ReceiVed: June 1, 2009; ReVised Manuscript ReceiVed: NoVember 17, 2009
QM/MM methods are widely used for studies of reaction mechanisms in water and protein environments. Recently, we have developed the VB/MM method in which the QM part is implemented by the ab initio valence bond (VB) method. Here, we report on further improvement of the VB/MM method which makes it possible to use the method for reactivity studies in systems where the QM and MM parts are connected by covalent bonds followed by first ab initio VB study of reactivity in proteins. We implemented a simple link atom scheme to treat the boundary interactions. We tested the performance of the link atom treatment in combination with the VB/MM method on an SN2 reaction and found it to be sufficiently accurate. We then used the VB/MM method to study the SN2 reaction in haloalkane dehalogenase (DhlA). We show that the predicted reaction barrier heights are in good agreement with estimated experimental values, thereby validating the method. Finally, we analyze the reaction energetics in terms of contributions of the VB configurations and conclude that such analysis is instrumental in pinpointing the essential features of the catalytic mechanism. Introduction A proper description of reactions in biomolecular systems requires the use of quantum mechanics/molecular mechanics (QM/MM) approaches. These methods enable the modeling of reactive biomolecular systems with satisfying accuracy and reasonable computational effort. Molecular orbital (MO) and semiempirical approaches dominate the QM part in most QM/ MM methods.1-21 However, differing from quantum studies of small reactive systems, there are also many QM/MM studies which are based on the valence bond (VB) approach and a variety of methods utilizing VB concepts have been developed.14,22-35 One such method is the empirical valence bond (EVB).14,22-24 In this method, the energies of the diabatic (resonance) states are represented by a classical potential which incorporates interactions of the charges with the environment. The EVB energies are then calibrated to reproduce experimental or ab initio data.14,22-24 The EVB benefits from the abilities of VB to describe chemical reactivity and was applied to various systems (e.g., refs 2 and 36-38). The success of the method in modeling reactions and its contribution to the understanding of reactivity led to the development of several extensions, trying to improve the description of the global shape of the potential energy surface (PES)25-31 or trying to generalize it to calculations of various VB states.32-35 A different VB based methodology that combines QM with MM is the molecular mechanics valence bond (MMVB) method39-41 where some active orbitals define the QM part of the system while the rest is described by MM. This method, however, is currently restricted to hydrocarbons. The weakness † Originally submitted for the “Walter Thiel Festschrift”, published as the October 29, 2009, issue of J. Phys. Chem. A (Vol. 113, No. 43). * Corresponding author. E-mail:
[email protected]. Phone: +9722-675-8696. Fax: +972-2-675-7076. ‡ Affiliated with the David R. Bloom Center for Pharmacy at the Hebrew University.
of this method and the methods described above is their dependence on parametrization of the QM part. Other methods utilizing VB involve the molecular orbital valence bond (MOVB) method42-44 which is based on the blocklocalized wave function approach to define the diabatic states45,46 and its extension to the effective Hamiltonian MOVB (EHMOVB) where a calibration process from EVB is adopted.44 Additional methods include the constraint or frozen density functional (CDFT and FDFT) where constraints are used to localize the electron density on various prechosen fragments, thereby describing diabatic states.47-50 Finally, two ab initio VB based methods were developed to describe reactions in solutions: VBPCM which combines VB with the polarizable continuum model51 and VBSM which combines VB with parameters determined for the SM6 solvation model.52 All of these unparameterized VB based methods, however, are either not suitable to describe reactions in biological systems or were never tested on such systems. Recently, an ab initio VB/MM method has been developed by us.53-56 The influence of the environment on the diabatic states, in this method, is considered by utilizing a point-charge (mechanical) embedding scheme. Mixing of these states then results in a proper polarization of the wave function.53,55 The present paper describes a further development of the ab initio VB/MM method in which treatment of bonded QM-MM interactions was introduced, thus allowing calculation of reactions in proteins. Bonded QM/MM interactions occur when a molecule is divided between the QM and MM regions, resulting in covalent bonds across the QM/MM interface. Models that treat these dangling bonds can be classified into three types. The first type concerns the link atom models, where a link atom is used in the QM calculation to cap the dangling bonds and saturate their free valency (see, e.g., refs 3, 15, and 57-63). The second type involves the connection (or boundary) atom models, where the boundary atoms interact with the other QM atoms as specially parametrized QM atoms and with the MM atoms as standard
10.1021/jp905143d 2010 American Chemical Society Published on Web 01/22/2010
VB/MM Protein Landscapes
J. Phys. Chem. B, Vol. 114, No. 6, 2010 2213
carbon MM atoms.64-66 The third type includes the hybrid (or localized) orbital model where hybrid orbitals are used to cap the QM region and are then kept frozen throughout the QM calculation.1,67-69 Both the hybrid orbital and the connection atom schemes are somewhat more rigorous than the link atom schemes. Nevertheless, various studies have shown that the results obtained by the link atom models are generally comparable with those obtained by both hybrid orbital models70-72 and the connection atom models,73 suggesting that the simple link atom schemes are sufficient. There are several variants of link atom schemes, differing from each other mainly in the treatment of the link atom positioning and the charge distribution on both the link atom and its neighboring MM atoms (e.g., refs 3, 15, and 57-63). Various studies have examined the effect of different link atom schemes.61-63,71-73 It was shown that the different schemes yield very similar structures and that the deviations between the relative energies and in particular reaction barriers are small on an absolute scale.61-63,71-73 This implies that the results are not very sensitive to the exact treatment of the QM/MM boundary, possibly due to error cancelation.62,63 Here, we present results of an ab initio VB/MM study of the SN2 reaction in the haloalkane dehalogenase protein in which bonded QM/MM interactions are treated by a simple link atom scheme. To the best of our knowledge, this is the first ab initio VB study of a reaction in protein environment. The paper is organized as follows. The first part of the paper describes the method with the assumptions made and their justification. The performance of the link atom based frontier treatment is then evaluated. Finally, the method is used to simulate the nucleophilic substitution reaction in haloalkane dehalogenase (DhlA) as a test case. It is shown that the method gives quantitatively correct predictions and provides further insight into the origin of enzyme catalysis. Methodology VB/MM Methodology Description. Details of the VB/MM method are presented elsewhere.53,55 Thus, here we describe the method only briefly with emphasis on the new developments. In VB theory, the wave function is expressed as a linear combination of multiple relevant VB configurations Φi74 N
ΨVB )
∑ ciΦi
(1)
i
The potential energy in VB/MM is, therefore, the lowest solution of the following equations:
HC - εSC ) 0
(2)
where H and S are the N × N diabatic Hamiltonian and overlap matrices, respectively, spanned by the N diabatic VB configurations describing the system and C is the N-dimensional coefficients vector. VB/MM like any other QM/MM method treats a small region where chemical change occurs with ab initio VB and includes the influence of the surroundings at the MM level. The Hamiltonian is therefore given by
HVB/MM ) HVB(I) + HMM(O) + HVB-MM(I, O)
(3)
where HVB(I) is a purely quantum VB Hamiltonian of all of the atoms in the inner region, HMM(O) accounts for the energy of
all of the atoms in the outer region determined classically by empirical force field, and HVB-MM(I,O) is the VB-MM coupling part of the Hamiltonian. HVB-MM(I,O) is responsible for all of the interactions between the quantum and classical atoms. Its evaluation, due to the multiple state character of the VB wave function, is not trivial and involves several steps and assumptions summarized here briefly. The interactions of the VB atoms with the environment are calculated separately for each VB configuration. Since these interactions differ for the different VB configurations, the interaction of the coupling term between the various VB configurations with the environment is not straightforward and is evaluated as the average interaction between each two VB configurations multiplied by their overlap as follows:
1 HVB-MM,ij ) (HVB-MM,ii + HVB-MM,jj)Sij 2
(4)
where HVB-MM,kk are the coupling term of VB configuration Φk with the environment. We note with this respect that eq 4 is based on the assumption that both the overlap integral, Sij, and the reduced resonance integral, βij, between the VB configurations are not affected by the environment. Similar to any QM/MM method, each of these coupling terms, HVB-MM,kk, includes bonded, VdW, and electrostatic interactions. b VdW elec HVB-MM,kk(I, O) ) HVB-MM,kk + HVB-MM,kk + HVB-MM,kk (5)
VdW interactions between the VB atoms and the environment are calculated classically. Similarly, electrostatic interactions are treated at the point charge (mechanical) embedding level for each individual VB configuration. That is, the wave function of each VB configuration in the QM region is translated into atomic point charges and the interaction with the MM region is then calculated using classical equations. We note with this regard that, as demonstrated earlier, mixing of the different VB configurations that follows guarantees polarization of the system’s wave function due to the environment (for more details, see refs 53 and 55). Thus far, the method was not designed to handle systems with bonding interactions at the VB-MM interface and the overall interaction with the environment was represented only by VdW and electrostatic interactions. Here, we introduce the ability to account also for covalent interactions between the VB and MM atoms. Link Atom Description. Previous studies have shown that the choice of link atom scheme is generally as good as other schemes.70-73 Furthermore, it was shown that different variations of the link atom scheme result with comparable relative energies.61-63,71-73 Since our interest is mainly studying reaction profiles, the simplest scheme seemed sufficient and was therefore employed. More specifically, a hydrogen atom is introduced along the broken VB-MM bond at an appropriate fixed distance from the QM atom.61,70 Two kinds of electrostatic interactions are excluded in order to guarantee balanced interactions. The first involves electrostatic interactions between MM atoms and the link atom, as the link atom is artificial and inclusion of its interactions with the environment is unjustified. The second involves electrostatic interactions between the VB atoms and MM atoms which originally are involved in bonding interactions (bonds, angles, and dihedrals) with any of the VB atoms (hereafter MMb) or MM atoms that complement the MMb group
2214
J. Phys. Chem. B, Vol. 114, No. 6, 2010
of atoms to a neutral group. Otherwise, all electrostatic interactions between the partial charges of the VB atoms and the MM atoms are included. The success of the scheme primarily depends on the manner in which the system is divided into VB and MM regions. The QM part is chosen in such a way that it carries a zero or whole charge while the remaining MM part is such that the partial charges of atoms in the MMb group are small in absolute value, leading to a total sum which is close to neutral. The bond being cut should be nonpolar and should not be involved in conjugative interactions. We emphasize that this scheme might not be the best scheme for certain specific properties (e.g., proton affinities, deprotonation energies),62,63 for which case other schemes can be implemented. Computational Details Quantum VB calculations were performed with the Xiamen VB (XMVB) program package,75,76 while one-electron integrals were obtained by the Gaussian 03 program package.77 Classical calculations including dynamics and all force-field related calculations utilized a modified version of the MOLARIS package.78,79 An interface program that communicates between the two programs was used. The VB/MM method describes a reaction by using several valence bond (VB) configurations. The method allows these VB states to mix and change their relative contribution along the reaction coordinate. In this study, the SN2 reaction of methyl chloride with aspartate in haloalkane dehalogenase was modeled in the gas phase, in solution, and in the protein. The system involves four electrons over three centers (4e/3c) and is defined therefore by six VB configurations, 1-6. Only four VB configurations were considered, reactant and product covalent R P and φcov , respectively (drawings 1 and 2), configurations, φcov the tri-ionic (ionic hereafter) configuration, φion (drawing 3), and the long bond configuration, φLB (drawing 4). Configurations 5 and 6, where either chlorine or the carboxylate group are positively charged at the presence of a carbanion, are expected to be very high in energy and contribute negligibly, and were thus eliminated from the calculation. VB gas phase calculations involved the VB self-consistent field (VBSCF) approach where the orbitals shared by the various configurations are optimized.80 These calculations were then used to get both the solution or protein results at the VBSCF/ MM level. The inner core electrons and the out-of-plane (πtype) electrons were frozen at the Hartree-Fock level. The remaining 24 and 32 valence electrons for the link atom and full model calculations, respectively, were explicitly included in the VB calculations. VB orbitals were all strictly localized on the respective fragments (RCO2, CH3, and Cl).
Sharir-Ivry et al. geometries were optimized at the MP2 level of calculation and the system was restrained to a CV symmetry. In the studies of the catalytic effect, additional restraints on the distance between the attacking oxygen and the leaving chloride and on the coordinates of the carboxylate group were placed to simplify the calculation in the protein. The remaining coordinates in this case were optimized at the HF level. We note in this respect that all of these restraints were shown to have a minor effect both structurally and energetically (see Table 1S in the Supporting Information). The link atom was fixed at a distance of 1.1 Å from the quantum boundary atom along the connecting bond. This distance which was obtained by a restraint optimization represents an average distance that would be obtained if the scaled bond distance approach would be employed.81 The classical simulations presented here were done using the ENZYMIX module of the program package MOLARIS.78,79 Spherical systems of radius 18 Å were used with the boundaries restrained by the SCAAS model.82 These surface restraints contribute to a proper behavior of the system’s density, calculated to be 0.97 gr/cm3 for the aqueous solution. Long range effects were treated by the local reaction field (LRF) long-range treatment.83 The starting coordinates used in the simulations of DhlA correspond to the 2DHC code taken from the Brookhaven Protein Data Bank (PDB)84 where dichloroethane was replaced by methylchloride. The initial distance between the attacking oxygen and the leaving chlorine atom was set to the corresponding values in the crystal structure. The calculations utilized the ENZYMIX force field.78 A detailed description of the force-field parameters for the reacting fragments at the various VB configurations is given in Table 2S and Schemes 1S and 2S in the Supporting Information. The reaction profile was calculated by following the potential of mean force PMF protocol,85-88 fixing the coordinates of the reacting fragments while relaxing the rest of the system. Each point along the reaction path was first relaxed for 50 ps with fixed coordinates for the reacting fragments and with some initial guess for the wave function (taken, for example, from a gasphase VB calculation or from results obtained for an adjacent geometry). A single point calculation using VBSCF/MM was then carried out and the resulting new wave function introduced as a guess for a new relaxation of the solvent/protein. The process was repeated until the changes in the weights of the various VB configurations were below 0.05. Once converged, each point along the reaction path was simulated for 80 ps in solution and for 50 ps in the protein. Geometry configurations were collected every 1 ps, and their respective energy was calculated at the VBSCF/MM level. All of the simulations were carried out at 300 K. Finally, the results were analyzed using the PMF procedure combined with the free energy perturbation (FEP) theory,89,90 to obtain the free energy profile of the reaction. Solvation energies of the various diabatic states as well as their contribution to the solvation of the overall adiabatic state both in solution and in the protein were determined by exploiting the linear response approximation (LRA) approach.91,92 The LRA treatment provides a good estimate for the free energy associated with the change between two potential energy surfaces U1 and U2 by eq 6:
1 ∆G(U1 f U2) ) (〈U2 - U1〉1 + 〈U2 - U1〉2) 2 The quantum calculations employed the 6-31G(d) basis set. To study the effect of the link atom model on the results, the
(6)
where 〈 〉j designates an average over trajectories propagated on Uj. Solvation free energy involves the interaction of the
VB/MM Protein Landscapes
J. Phys. Chem. B, Vol. 114, No. 6, 2010 2215
reacting fragments with their environment (water/protein). Therefore, the two potential energy surfaces, U1 and U2 in that case, represent the energy of the whole system (reacting fragments and their environment) when the partial charges of the reacting fragments are considered, U1, and when they are set to zero, U2. Hence, in order to get the solvation energy, we have to run two simulations, one for each electronic state. Using LRA enables decomposing the free energies to their individual additive contributions, allowing analysis of the role of the various diabatic states. The LRA calculations involved averaging of 80/50 different structural configurations of the water/protein, respectively, collected for each geometry of the reacting fragments. This was done by continuously equilibrating the system at each solute geometry and performing a VB/MM calculation after every 1 ps in simulation that spanned an overall 80/50 ps. All of the simulations were done at 300 K with a 1 fs time step.
TABLE 1: Reaction Barriers (kcal/mol) for the SN2 Reaction Described in eq 7 Both in Gas Phase and in Aqueous Solution (Comparison of the Link Atom Scheme with Full Calculation)
Results
R φcov P φcov
Our goal is to study enzyme catalysis with the ab initio VB/ MM methodology. The enzyme we chose is the widely studied haloalkane dehalogenase (DhlA).93-101 DhlA catalyzes the hydrolytic cleavage of carbon-chlorine bonds in a wide range of haloalkanes, turning them into their corresponding alcohols. The reaction chosen for this study is the first step in the DhlA reaction which involves a nucleophilic substitution of the haloalkane by the carboxylate group of the Asp-124 residue followed by formation of a covalent bond with the alkane and resulting in an alkylated enzyme intermediate and a halide ion. Validity of Link Atom Scheme. Before embarking on the study of the reaction in the protein environment, it is necessary to examine the effect of the link atom scheme on VB/MM results. This can be easily accomplished by carrying out calculations on a smaller system which is amenable to treatment by both VB and VB/MM methods. Such a system will facilitate direct comparison of results that were obtained, on one hand, by treating the whole system at the quantum mechanical level (termed full hereafter), and, on the other hand, by dividing the system into MM and QM regions connected by a covalent bond (referred to as link atom). Equation 7 shows the simplified system that was used for the evaluation of the link atom scheme. It is a model of the first step in the conversion of haloalkanes into alcohols by haloalkane dehalogenase. The catalytic residue Asp-124 is represented by propanoate that contains the carboxyl group, and we chose methyl chloride as a representative haloalkane substrate:
CH3CH2COO- + CH3Cl f CH3CH2COOCH3 + Cl(7) The boundary between the QM and MM regions was chosen to be the CC bond that connects the ethyl and carboxyl groups. Table 1 summarizes the comparison of reaction barriers obtained by the two methods, full and link atom, both in the gas phase and in aqueous solution. Looking at the results, we find that the reaction barriers are essentially independent of the boundary treatment both in gas phase and in aqueous solution. This is encouraging, since it means that the accuracy of the VB calculations did not deteriorate upon introduction of a link atom. Furthermore, the difference in the barrier height prediction between the full and link atom treatments in the gas phase (0.4 kcal/mol) is very similar to that obtained in aqueous solution
boundary treatment
∆gqg
∆gwq a
full linkatom
11.9b 11.5a
18.2 18.7
a Free energies are calculated at the [VBSCF/6-31G(d)//MP2/ 6-31G(d)]/MM level using PMF. b Energies are calculated at the VBSCF/6-31G(d)//MP2/6-31G(d) level.
TABLE 2: Weight of the Various VB Configurations Both in Gas Phase and Solution (Comparison of Link Atom Scheme with Full Calculation)a gas full
φion φLB
water link atom
full
link atom
RS
TS
RS
TS
RS
TS
RS
TS
0.615 0.005 0.378 0.002
0.290 0.137 0.554 0.020
0.612 0.007 0.379 0.002
0.291 0.134 0.555 0.019
0.646 0.003 0.350 0.001
0.240 0.134 0.612 0.015
0.648 0.002 0.349 0.001
0.261 0.113 0.612 0.013
a RS and TS geometries are taken from the VB profile at each environment.
(-0.5 kcal/mol). Such a small difference can be expected due to different sampling in the gas phase and aqueous solution. Moreover, the small magnitude of the difference suggests that errors associated with link atom treatment are insensitive to polarization by solvent. An important and unique property of the VB methodology is the wave function of the system under scrutiny. Thus, it is important to examine the effect of the model on the wave function. The comparison of the weights of the VB configurations together with the reaction energetics provides a more stringent test of the method’s accuracy than reaction energetics alone. Table 2 compares the weights obtained by the two treatments for the VB configurations in the ground state and transition state geometries, both in gas phase and in aqueous solution. Clearly, the differences between the two schemes are very small in gas phase calculations. In solution, we observe somewhat larger differences but still sufficiently small lest to warrant further investigation. The higher magnitude of weight discrepancy in solution is most likely caused by issues related to convergence of sampling during the simulation. We note here that, taken as a whole, the weights in gas phase are similar to the corresponding weights obtained in aqueous solution. This however is in agreement with the effect of the aqueous solution on the weights on a similar system, as demonstrated by Shaik and co-workers.102 Overall, the inclusion of a link atom in the VB/MM protocol did not cause any major errors in this particular system and we therefore feel confident that this methodology is suitable for the studies of haloalkane substitution in the protein environment. The Catalytic Effect. The reaction mechanism of DhlA has been widely studied by many different approaches; however, there are no ab initio VB studies. Here, we fill this gap and present the results of a VB/MM study for the first step of the reaction. The energy profiles along the reaction coordinate were calculated in gas phase, water, and protein environments with methyl chloride (MC) as the substrate. The reference model system of eq 7 was used in gas phase and water. Table 3
2216
J. Phys. Chem. B, Vol. 114, No. 6, 2010
Sharir-Ivry et al.
TABLE 3: Barrier Heights (in kcal/mol) for the Methylchloride SN2 Reaction in Gas Phase, Aqueous Solution, and in the DhlA Enzyme along with Differences in Barrier Height due to Change of Environment Expa VBSCFe
∆Eqg
∆gwq
∆gqp
q ∆∆gwfg
q ∆∆gpfw
