Insights on the origin of catalysis on glycine N-methyltransferase from computational modelling

The origin of enzyme catalysis remains a question of debate despite much intense study. We report a QM/MM theoretical study of the SN2 methyl transfer reaction catalyzed by a glycine N-methyltransferase (GNMT) and three mutants to test whether recent experimental observations of rate-constant reductions and variations in inverse secondary α-3H kinetic isotope effects (KIEs) should be attributed to changes in the methyl donor-acceptor distance (DAD): Is catalysis due to a compression effect? Semiempirical (AM1) and DFT (M06-2X) methods were used to describe the QM subset of atoms, while OPLS-AA and TIP3P classical force fields were used for the protein and water molecules, respectively. The computed activation free energies and KIEs are in good agreement with experimental data, but the mutations do not meaningfully affect the DAD: Compression cannot explain the experimental variations on KIEs. On the contrary, electrostatic properties in the active site correlate with the catalytic activity of wild type and mutants. The plasticity of the enzyme moderates the effects of the mutations, explaining the rather small degree of variation in KIEs and reactivities.


INTRODUCTION
Glycine N-methyltransferase (GNMT) is an S-adenosyl-L-methionine (SAM)-dependent enzyme that catalyzes the transformation of glycine into sarcosine (see Figure 1a). 1 SAMdependent methyltransferases are inhibited by S-adenosyl-L-homocysteine (SAH), and it is thought that GNMT can be a regulator of the cellular SAM/SAH ratio, playing a key role in other methyl transfer reactions. 2 Thus, the activity of GNMT has been related with the concentrations of SAH and homocysteine in plasma which, in turn, has been linked to dementia and to Alzheimer's disease. 3 Nevertheless, although the catalytic mechanism of GNMT is well known based on kinetic and crystallographic studies, 4,5,6 the origin of the rate enhancement in the chemical step, a SN2 methyl transfer, is still a question of debate. The most accepted hypothesis of the origin of the rate enhancement in this (or any other) enzyme is the relative stabilization of the transition state (TS) with respect to the reactant state (RS), as originally proposed by Pauling. 7 In other words, the reaction catalyzed by an enzyme has a lower freeenergy barrier than the counterpart reaction in aqueous solution. Based on molecular simulations, this reduction of the free-energy barrier has been attributed to electrostatic preorganization of the protein that stabilizes the TS with respect to the RS. 8 Nevertheless, for enzyme catalyzed SN2 reactions, it has been suggested that specific protein fluctuations might reduce the donor-acceptor distance (DAD), thus diminishing the potential-energy barrier height and/or width and enhancing the rate by increasing the number of reactive trajectories over and through the barrier. 9,10,11,12,13,14 Obviously, in reactions where the transferred moiety is not a light particle, as in the case of GNMT, a large contribution from tunneling is not expected. The idea of reduction in the DAD has its origins in Schowen's "compression hypothesis" 15 for enzymic methyl transfer. This was based on experimental observations that the secondary kinetic isotope effect (2° KIE) k(CH3)/k(CD3) for enzyme-catalyzed methyl transfer was much more inverse than for the uncatalyzed reaction in solution, suggesting a "tighter" TS for the enzymic reaction. Mechanical compression by the enzyme destabilizes the RS more than the TS, resulting in a lower free-energy barrier. Whether or not this idea is valid has since been argued based on theoretical 16,17,1819,20 and experimental studies. 21 Thus, the controversial contribution of protein motions to modulate the DAD, originally proposed in enzymatic hydride transfer reactions, 9,10,11,12,13 has been extended to GNMT 22 and other related methyltransferases such as catechol O-methyltransferase (COMT). 23,24 a. Chemical step b. Full enzymatic process Figure 1. a) Schematic representation of the chemical step of the GNMT catalyzed reaction. b) Full catalyzed mechanism for the transformation from SAM into SAH catalyzed by GNMT. The particular isotope effect is indicated on each step.
In order to quantify and rationalize the origin of catalysis in enzymes, calculations can be carried out in the framework of Transition State Theory (TST). 25,26 Within this context, and including corrections for possible tunneling contributions and dynamic effects, 27,28 computed rate constants can be directly compared with those from kinetic experiments. In particular, computational QM/MM methods that combine quantum mechanics (QM) and molecular mechanics (MM) potentials have emerged as the most appropriate strategy for study of chemical reactions in very large flexible systems such as enzymes.
Several studies based on hybrid QM/MM methods have been carried out in our laboratories on enzymatic methyl transfer reaction catalyzed by COMT and, according to our studies, the origin of catalysis should be attributed to the electrostatic stabilization of the TS, 17,29,30 with a low coupling between the protein motions and the reaction coordinate on the chemical step of COMT. 31 In addition, determination of 2º KIEs for COMT catalyzed reaction did not support compression effects. 17 Regarding GNMT, our previous study on the rat-liver wild type enzyme (WT) showed how the binding energy is maximal for the TS 32 and the averaged DAD for the computed TS structures was actually smaller in solution than in the enzyme. These results rule out the compression along the donor-acceptor axis for these two methyl transferases. However, these conclusions have subsequently been questioned. The experimentally measured reduction of the catalytic efficiency of GNMT when active-site residue Tyr21 was replaced by a series of mutants, as well as the change in the 2º -3 H KIEs, was interpreted by Zhang and Klinman as a catalytically relevant reduction in the methyl DAD. 22 Here we study the methyl transfer reaction catalyzed by WT rat-liver GNMT and its Y21F, Y21A and Y21G mutants by means of variational TST, and compare with the counterpart reaction in aqueous solution in order to shed light on the source of catalysis.

COMPUTATIONAL METHODS.
Set up of the models. The molecular system was prepared based on a crystal structure of WT rat-liver GNMT (PDB ID 1NBH) at 2.00Å resolution. 5 Since the standard pKa values of ionizable groups can be shifted by local protein environments, 33 an accurate assignment of the protonation states of all these residues at pH = 7 was carried out. Recalculation of the pKa values of the titratable amino acids has been done using the empirical PROPKA 2.0 program of Jensen et al. 34 Accordingly, Asp62 was protonated in OD2 position, His142, His264 and His286 were protonated in δ position, His58, His153, His174, His214, His245 were protonated in ε position, and finally His55 was protonated in both δ and ε position respectively.
Subsequently, in order to neutralize the system, one counterion (Cl -) was placed into the optimal electrostatic position around the protein (that is where the potential reaches its maximum positive value). Afterwards, a series of optimization algorithms (steepest descent, conjugate gradient, and L-BFGS-B) 35 were applied. To avoid a denaturation or loss of structure prior to solvation of the protein, all the heavy atoms of the protein were restrained by means of a Cartesian harmonic umbrella with a force constant of 1000 kJ·mol -1 ·Å -2 . The protein was placed in a box of pre-equilibrated water molecules (100  80  80 Å 3 ), using the principal axis of the protein-inhibitor complex as the geometrical center. Any water with an oxygen atom lying in a radius of 2.8 Å from a heavy atom of the protein was deleted. The geometries of the remaining water molecules were then optimized. Later, an initial 500 ps classical MD simulation (at temperature 293 K), carried out to relax the system without significant changes in its geometry, was performed using the AMBER 36 and TIP3P 37 force fields for the protein and water molecules, respectively, as implemented in the fDYNAMO library. 38,39 The same force fields were used to describe the MM region in the subsequent hybrid QM/MM MD simulations. In this case, a small part of the system, consisting of the full SAM and glycine as depicted in Figure 2, was described by QM using either the AM1 40 semiempirical Hamiltonian or the M06-2X hybrid functional 41,42 with the standard 6-31+G(d,p) basis set. Due to the large number of degrees of freedom, any water molecule 20 Å away from any of the atoms of the substrate and cofactor was kept frozen in the remaining calculations. Cutoffs for the non-bonding interactions were applied using a force-switching scheme, within a range radius from 14.5 to 16 Å.
The same computational protocol was followed to prepare systems corresponding to the Y21F, Y21G and Y21A mutants. In these cases, the replacement of the original tyrosine residue on position 21 by Phe, Gly and Ala was performed on the X-ray structure and, consequently, the systems were equilibrated by hybrid QM/MM optimizations and MD simulations. Finally, a reference system in aqueous solution was prepared. The full SAM and substrate Gly molecules were embedded in a box of water molecules (100 × 80 × 80 Å 3 , 19818 solvent water molecules in total) treated also by the TIP3P force field, as implemented in the fDYNAMO library. The same non-binding interaction conditions were applied to prepare the system in solution. Water molecules were equilibrated by means of hybrid AM1/MM MD in the reactant complex.
Long MD simulations of 10 ns were finally done for each of GNMT variants under NVT conditions with T = 293 K. In order to run long MD simulations NAMD software was used. 43 All models were described at MM level of theory using AMBER force field for enzyme and Gly substrate and contraions, TIP3P force field for water molecules. Missing parameters of SAM, determined using antechamber module available in AmberTools, 44,45 are reported in the Table S1 of Supporting Information. Because compression is a key aspect of the present problem, it is interesting to wonder whether the results are sensitive to employing an isobaric NPT ensemble rather than the NVT ensemble. Simulations at higher pressure have shown that donor-acceptor distances in the hydride transfer reaction catalyzed by morphinone reductase can be reduced by 0.1-0.2 Å when the pressure is increased from 1 bar to 2000 bar. 46 However, the impact on the free energy barrier is expected to be small because the force constant for donor-acceptor stretching in the RS is quite small. In the case of GNMT a similar trend could be expected. Instead, the effect of a similar pressure increase on the TS geometry must be quite limited considering that the force constant for donor-acceptor stretching is larger in the TS. In fact, simulations of the reactant ternary complex of the WT enzyme under NPT conditions, at 1 and 1000 bar, rendered average structures with DADs statistically indistinguishable from the values derived from the NVT MD simulations (see Table S2 in Supporting Information).
Free energy surfaces. In order to get the free-energy landscape of the catalyzed reaction, potentials of mean force (PMFs) were traced along the selected reaction coordinate using the weighted histogram analysis method (WHAM) 47 combined with the umbrella sampling approach, 48 as implemented in fDYNAMO. In the present study, two-dimensional (2D) PMFs were generated by explicitly controlling the antisymmetric combination of distances defining the methyl transfer from the sulfur atom of SAM to the nitrogen atom of Gly, d(S-C) -d(C-N), ξ1, and the distance between the donor and acceptor atoms, S and C, as a second reaction coordinate, ξ2. The values of the variables sampled during the simulations are then pieced together to construct a distribution function from which the PMF is obtained as a function of the distinguished reaction coordinate (W(ξ)). Then, the activation free energy of a chemical step can be expressed as: 49 where the superscripts indicate the value of the reaction coordinate at the reactants (R), and at the TS ( ‡), and Gξ (ξ R ) is the free energy associated with setting the reaction coordinate to a specific value at the reactant state. Then, the quasiclassical activation free energy , as appearing in Eq. 1, is calculated along the reaction coordinate , 50,51 as the sum of ∆ ‡ ( ) and the correction term due to the quantized nature of molecular vibrations (mainly zero-point energies). 50,52 Because the generation of PMFs requires generation of a very large number of structures from QM/MM MD calculations, inevitably, we are restricted to the use of a semiempirical Hamiltonian (AM1) in this work. However, in order to improve the quality of our free-energy surfaces, based on the work of Truhlar and co-workers, 53,54,55 a spline under tension 56,57 is used to interpolate a correction term at any value of the reaction coordinates ξ1 (and of ξ2 in the case of 2D PMFs) selected to generate the free-energy surfaces. 58,59,60 In this way a new continuous energy function is constructed in order to obtain corrected PMFs: where S refers to the spline function, and its argument is a correction term evaluated from the single-point energy difference between a high-level (HL) and a low-level (LL) calculation of the QM subsystem. Herein the the M06-2X functional with the 6-31+G(d,p) basis set, as suggested by Truhlar and co-workers, 41 Hessians have been subjected to a projection procedure to eliminate translational and rotational components, which give rise to small nonzero frequencies, as previously described. 64 Thus, it has been assumed that the 3N − 6 vibrational degrees of freedom are separable from the 6 translational and rotational degrees of freedom of the substrate.
Finally, as mentioned in the introduction, a qualitative change in KIEs is not expected if quantum tunneling effects were to be incorporated into the calculations, taking into account the mass of the transferring particle (a methyl group). A test was done with the Wigner method confirming this prediction: an average nuclear quantum correction of 1.003 was obtained for both isotopologs, cancelling out the effect on the KIE.

RESULTS AND DISCUSSION
The free-energy surfaces (FESs) generated from AM1/MM simulations for the reaction in water, WT GNMT and its Y21F, Y21A and Y21G mutants are shown in Figure 3 5 respectively. The increases for Y21A and Y21G are higher (4.2 and 5.9 kJ·mol -1 relative to WT) but less than the experimental differences (14 and 13 kJ·mol -(see Figure 4) or at AM1/MM show similar DAD values and the same conclusion is derived from the analysis of RS DADs (see Table S3 in Supporting Information).  Table S3

Isotope Effects.
Our results show a good qualitative agreement between the experimental and computational 1 H3/ 3 H3 binding isotope effects (BIEs) for all four GNMT variants (see Table 1). These normal BIEs (values larger than unity) derive from stiffer vibrational modes associated with these H atoms in aqueous solution than in the enzyme active site. Geometrical analysis of the optimized structures shows that these H atoms make more and tighter interactions with water molecules in solution than with the residues in the GNMT:SAM binary complex of any of the protein variants (see Figures S4 and S5 in Supporting Information). This is reasonable since the active site must present a cavity before it can accommodate the binding of glycine. Table 1. 1º 12 C/ 14 C KIEs and 2º 1 H3/ 3 H3 KIEs between the enzyme:SAM binary complex and the TS ( theor KIE b in Figure 1b). 12  The calculations do indeed predict inverse 2º 1 H3/ 3 H3 KIEs for all the proteins and slightly less inverse in aqueous solution as measured experimentally, but no meaningful trend can be observed (within the estimated 1 uncertainty) between the GNMT variants. The origin of the inverse 2º 1 H3/ 3 H3 KIEs must be related with an increase of the force constants associated to the methyl group hydrogen atoms from the GNMT:SAM binary complex to the TS (see Table   S7 in Supporting Information), which is also reflected in the evolution of the C-H bond lengths diminishing from the GNMT:SAM:Gly ternary complex (RS) to the TS (see Figure S7 in Supporting Information). Among the four proteins, the WT stabilizes the TS better than the mutants (as reflected by the predicted and experimentally deduced free energy barriers) which can be related a priori with the loss of the stabilizing interaction between the methyl group and the hydroxyl group of Tyr21 in the WT after mutation. However, as shown in Figure 4, the interaction between a H atom of the methyl group and His124 appears to be stronger in mutants (shorter inter-atomic distance in the three mutants). In addition, a second H atom of the transferring methyl group interacts with the backbone carbonyl group of Gly137 in all the optimized TS structures with an O-C distance in general slightly shorter than in the GNMT:SAM binary complex (see Figure S5 in Supporting Information). Apparently these interactions are thus compensating for the loss of the interaction with Tyr21 and explaining why the 1 H3/ 3 H3 BIEs and 2º KIEs do not differ sensibly between WT and the mutants, within the standard deviations.
The computed 1º 12 C/ 14 C KIEs are normal and in an almost quantitative agreement with the experimental results, and the 12 C/ 14 C substitution produces no BIEs; but no experimental data were reported for this substitution to be compared with.
Electrostatic effects. The charges on the key atoms for donor, acceptor and the transferring methyl group (computed with the CHELPG method) 66 change in the same way along the reaction coordinate (r.c.) in aqueous solution and in WT GNMT ( Figure 5). The charge on the donor S atom varies from slightly positive in the RS to slightly negative in the TS, in both media. The acceptor N atom is always more negatively charged in the enzyme than in solution, and the sum of the methyl-group atomic charges is always more positive in the enzyme than in solution. So it seems that the protein environment favours a TS with larger separation of charges between the transferring methyl group and the accepting N atom. The evolution of the charges on the three mutants along the r.c., deposited in the Supporting Information (see Figure S8 in the Supporting Information), are almost coincident with the WT shown in Figure 5. It is interesting that the methyl-group charge does not reach its maximum in the TS but after passing this structure. This can be related not only to the evolution of the DAD discussed above, but to the evolution of the electrostatic properties of the different environments. Thus, as shown in Figure 6, the averaged electrostatic potential created by each of the different environments at the three key atoms changes from the RS (r.c. = -1.5 Å) to the point where DAD reaches its minimum (r.c. +0.5 Å). All environments appear to be equilibrated from this point on.
Consider the three key atoms. At the N atom the potential due to each environment changes from positive in the RS to negative in the PS. Since this atom is itself negatively charged in the RS and the magnitude of its charge decreases as the reaction proceeds, this means that all environments, from the point of view of this atom, stabilize the RS rather than the TS or PS.
Similarly, since the charge on the S atom evolves from slightly positive in RS to slightly negative in PS, and the potential becomes less negative as the reaction proceeds, all environments also stabilize the RS. However, the negative potential created by each environment at the C atom stabilizes the positive charge developed on the transferring methyl group in the TS. Interestingly, there is a linear correlation between this potential and the calculated free energy barriers (more negative potentials correspond to smaller barriers) which corroborates the electrostatic origin of catalysis in GNMT (see Figure 6F). The WT enzyme is the environment that generates the most negative potential at the TS thus stabilizing better the positive charge developed on the methyl group. This means that WT provides the most favourable environment within which to transfer the methyl group from S to N (see Table S8 in the Supporting Information for the full list of averaged electrostatic potentials). On the other hand, analysis of the contributions to the total electrostatic potential by residues supports the small differences observed in the four proteins. While the potential generated by residue 21 is less favourable to stabilize the positive charge developed in the methyl group in the TS after mutations, the potential generated by His142 is more favorable in the mutants than in the WT.
Comparison of the plots presented in Figure 6 reveals an important point. Whereas the change from RS to PS in the electrostatic potential due to the environment is quite different in water as compared to any of the proteins, these changes are all very similar for the WT enzyme and the three mutants. This accords with a higher energy of reorganization from RS to TS in aqueous solution, and with only small differences in enzyme reactivity consequent upon mutation of Tyr21. Thus, a priori, the catalytic effect of GNMT could be partially attributed to the generation of a non-polar cavity: since the reaction catalyzed by GNMT implies annihilation of charges from reactants to products, a polar medium (like water) would preferentially stabilize the reactants and consequently higher barriers should be expected. Nevertheless, our calculations reveal that this simplified picture is not complete, since the enzyme is generating a potential that facilitates the positively charged methyl transfer and the stabilization of the TS.
Catalysis in GNMT is not only due to the absence of a strong reaction field (as the one felt in water by the reacting system) but also to the existence of a preorganized electric field favouring the methyl transfer. Indications of this catalytic effect have been previously detected in the methyl transfer reaction catalyzed by COMT, 30

CONCLUSIONS
The first important conclusion derived from our QM/MM study is that the activation free energy of the reaction in solution is significantly higher than in any of the enzymes, which reflects the catalytic role of the proteins. Although the differences are small, the WT enzyme shows the lowest activation free energy barrier while Y21A and Y21G mutants lead to the highest values.
These results are in very good agreement with the values that can be derived from the experimentally measured rate constants by Takata et al. 5 and by Zhang and Klinman. 22 The analysis of the average DAD in RS and TS does not provide any indication about the compression effect as a source of catalysis since similar values are obtained in all environments.
In fact, structural analysis of the active site shows that Tyr21 cannot push the S donor atom towards the N acceptor atom.
The more inverse 2º 1 H3/ 3 H3 KIEs observed for the WT GNMT than for the mutants has been interpreted as being due to a "compactness/tightness" of its active site as compared with the mutated variants, which have diminished ability to achieve a constrained DAD. 22  The detrimental effect of substitution of Tyr21 on the electrostatic potential is partially compensated by His142 that, by approaching to the methyl group, generates a higher potential.
The changes on the electrostatic potential exerted by the different proteins on the substrate reflect that they are not static but dynamic macromolecules that have to evolve structurally from the RS to the TS, although at a low energy cost. The plasticity of the WT protein is responsible for the differences in the KIEs and rate constants after mutations being as small as they are.

SUPPORTING INFORMATION
The Supporting Information is available free of charge on the ACS Publications website at DOI: