Methods and strategies

Methods:
The reaction has been modeled using a QM/MM potential as implemented in AMBER 5.0[52] (Roar-cp[235] module). Roar-cp module is the result of coupling together SANDER (the standard AMBER minimizer) with mopac 7.0. This module has been developed in the Merz laboratory during the past decade. A modified version of this program has been used.

The QM/MM scheme used here is displayed in equation 1.62 in page . Hydrogen link atoms have been employed to cap exposed valence sites due to bonds which cross the QM-MM boundary. The link atom model used in this work is the same that is used in the original reference [68]. It means that the hydrogen link atom does not feel any interaction from the MM part at all. It has been observed that the model is better reproduced when the link atoms can interact with all the MM atoms excluding the MM host[80]. Doing so we would have avoided the artificial polarization that we eventually had found. In addition, the link atom is left free as an additional degree of freedom. It is permitted to move during optimizations and its coordinates are saved as any other atom. Despite other approaches freeze at 1.0Å the distance of the link atom, we think that in the optimization process this approach would cause problems on the convergence criterion.

As in the previous work on vinylglycolate[212], the QM part of the system has been represented at the PM3 semiempirical molecular orbital level, but for the magnesium ion the new AM1 parameters developed by Hutter et al.[234] have been used. Several recent theoretical studies on biological systems have successfully used PM3[32,33] and AM1 [31] methods (see reference [72] for a review of the applications).

As for the MM part of the system and the QM/MM van der Waals and electrostatic interactions we have utilized the AMBER force field by Weiner et al[51]. The full QM/MM energy computation has been carried out taking into account all the interactions, that is, no cutoff was applied for the non-bonded interactions during the most of the simulations (with the exception of the molecular dynamics run).

For the sake of comparison with the case of vinylglycolate, in the present work we have chosen the same QM/MM partition used in the Garcia-Viloca et. al. study (see figure 2.3 page and table 2.2 page ). This model has a QM subsystem with zero total electronic net charge, and 80 or 88 atoms, link atoms included, for propargylglycolate or mandelate complexes, respectively. This means that a total of 3963 atoms will constitute the whole QM/MM model for the mandelate case.

In QM/MM scheme, when cutting a covalent bond, the MM part of the cut residue can have atomic charges which sum is not an integer. In order to preserve integral charge in the MM region, in each partially QM modeled residue, the MM charges of the lateral chain atoms have been changed[236]^3.7.

At the moment of writing this thesis the method used in this chapter (carried out more than three years ago) may appear too rough and rather antiquated. This is a sign that fortunately this field and the computational power is changing very rapidly. However, we believe that the chemical insight that these calculations give is still useful. Setup of the model:
The 2.0 Å resolution structure of the complex of Pseudomonas putida Mandelate Racemase enzyme inactivated with (R)-$\alpha$-phenylglycidate (Protein Data Bank code 1MNS[233]) has been taken as the starting point for the calculations. From this crystallographic structure we have taken the Cartesian coordinates of the 2700 atoms that belong to enzyme residues and we have discarded the coordinates corresponding to the 209 crystallographic waters, with the exception of the water directly bound to the metal Mg ion in the active site. Instead of keeping the crystallographic waters, we have solvated the active site with a sphere of TIP3P[237] water molecules of radius 20 Å, centered on the magnesium atom. This sphere includes 202 water molecules that are submitted to a soft harmonic potential to avoid their moving away from the active center.

We have substituted the (R)-$\alpha$-PGA inhibitor found in the X-ray structure by the different substrates studied here, that is, (S)-propargylglycolate and (S)-mandelate, superimposing as many atoms as possible. Since no hydrogen atoms were determined by the crystallographic technique, they have been added using the EDIT module of AMBER[51]. However, we have not represented explicitly all the hydrogen atoms because we have used AMBER united atom model for the enzymatic residues that belong to the MM part of the system.

The protein complex has been neutralized by placing four Na$^+$ in positions of large negative electrostatic potential and far enough from the active site.

A 15 Å sphere was defined around the active site magnesium atom, and only residues within this sphere as well as the water molecules were allowed to move during the simulations, that results in 1299 moving atoms for the mandelate substrate case.

We have also carried out QM/MM molecular dynamics simulations in the NVT ensemble to find new possible structures that could represent the reactant and the product complexes of propargylglycolate and mandelate with the enzyme. Starting from the minimized structures, a QM/MM molecular dynamics simulation was performed with a 1 fs time step to heat the systems from 0 to 300 K over an interval of 6 ps, with atom velocities assigned from a Gaussian distribution every 2 ps in 100 K increments. The systems were equilibrated for an additional 10 ps run. The temperature was maintained by coupling to two thermostat chains (one for the MM region and the other for the QM region) within the Nose­Hoover chain temperature scheme [174,176]. We have used three thermostats for each chain and the program has automatically set the mass values of each of them. In this case the non-bond cutoff distance was 10 Å for the MM atoms whereas we used no cutoff for the QM atoms. The non-bond pair list was updated every 25 steps and the SHAKE algorithm[178] was used to constrain bond distances that imply hydrogen atoms. The structures obtained after the molecular dynamics simulation turned out to be very similar to those corresponding to the minimized complexes. This ensures that both structures, that is, before and after the molecular dynamics, are a good starting point for the study of the reaction mechanisms. Procedure for the energy profiles:
Starting from the model structure of the reactive complex for each different substrate, the QM/MM potential energy was minimized until the RMS gradient fell below 0.001 kcal/(mol$\cdot$Å), by means of the L-BFGS method[153] (see section 1.3.7.1 in page for a description of the method). The resultant minima were taken as the reference structure that models the reactant of the S$\to$R reactions. From each reference structure we have calculated the energy profiles of different possible mechanisms of the enzymatic reactions.

We have minimized the QM/MM potential energy along a suitable reaction coordinate for each of the reactions, with the convergence criterion described above.When the step consists basically in a proton transfer the distance between the acceptor atom and the hydrogen that is being transferred is taken as the reaction coordinate. When the reaction is a conformational change (inhibition of propargylglycolate in section 2.3.5) the appropriate dihedral angle is scanned. In any case, a penalty harmonic potential is applied to this reaction coordinate: U = K(r - r$_{eq}$)$^2$, where r$_{eq}$ defines every point of the corresponding energy profile, and K = 10 000 kcal/(mol Å$^2$), that permits the scanning of the reaction coordinate in small increments of 0.1 Å. In this section, the structure of maximum classical potential energy at every step of the mechanism will be considered as the transition state of that step.

The steps studied here for propargylglycolate and mandelate substrates are those described in the previous work for vinylglycolate. For a full treatment of a chemical reaction, whether in a gas or condensed phase, the dynamics of the process should be included. However, studies of the potential energy surfaces[160,84,161] governing enzymatic reactions can in themselves be mechanistically revealing, and are an essential precursor to more extensive dynamic studies.

First, we will present the results corresponding to propargylglycolate. Secondly, the case of mandelate will be analyzed. All reactions have been studied from the (S)-enantiomer to the (R)-enantiomer. We recall that the relevant bond distances along with the energetics of the different intermediates are displayed in the tables of section 2.3.8 page .