Combining two bond distances

NHC and CHN:
Starting from the equilibrated S structure the reaction coordinate R$_{NHC}$, as defined in table 4.2, is scanned by the different windows in the S $\to$ R direction. After this simulation the final equilibrated R structure is taken to scan the R$_{CHN}$ reaction coordinate in the R $\to$ S direction. Unless a major sampling is needed the reference reaction coordinate $Rc_0$ used in the umbrella potential is increased by units of 0.2 Å.

A harmonic restraining potential is often used as biasing function. Sometimes, in order to obtain an efficient sampling an additional biasing potential is used[261]. A cubic spline function fitted to a number of points along the reaction coordinate is added to the harmonic restraining potential.

$$U_{tot} = U_{QM/MM} + k(Rc-Rc_0)^2 + U_{spline}$$ (5.1)

This strategy would obtain an enhanced exploration of the low probability zones without increasing the harmonic constant. The drawback is that this additional biasing function should be ideally equal to the negative of the unknown PMF ( $-\Delta G(Rc)$) that is actually what we want to calculate. Then the points to which the spline is fitted must be determined by trial and error until obtaining the expected sampling.

We have used this additional biasing in the PMF for R$_{NHC}$ and R$_{CHN}$ reaction coordinates. In this case the spline function for an appropriate sampling was found easily thanks to the energy profile obtained in the preceeding chapters by optimization of the structures. In the rest of PMF calculations such trial error effort was not worthwhile and the sampling is only improved by increasing the restraining harmonic constant. In the present case the harmonic constant is 20.0 kcal/(molÅ$^2$) while for the rest of reaction coordinates is increased to 30.0 kcal/(molÅ$^2$).

In figure 4.4 we show the histograms for the probability distribution $P(Rc)$ for the R$_{NHC}$ and R$_{CHN}$ reaction coordinates. $P(Rc)$ is converted to the probability density $\rho (Rc)$ by removing the biasing effect in the posterior statistical analysis (see introductory section 1.5.2).

Figure 4.4: Representation of the several probability distributions P(Rc) collected for the different windows. Units are given in number of configurations

It is considered that the sampling has converged when the bin is visited more than approximately 500 times^5.2.

In figure 4.5 the potential of mean force is shown for the R$_{NHC}$ (S$\to$R) and R$_{CHN}$ (R$\to$S).

Figure 4.5: PMF profile for R$_{NHC}$(left) and R$_{CHN}$(right)

The scanning of the reaction coordinate R$_{NHC}$ (on the left of figure 4.5) reaches a zone after the proton transfer from substrate to Lys166 (R $_{NHC}\sim 1.0$) where this coordinate is unable to describe the carbon configuration inversion. After this region the simulation explores a zone that does not belong to the reaction path anymore. At this point Lys166 moves away from the substrate, leaving the stereogenic carbon in S configuration.

The PMF calculation using R$_{CHN}$ (the profile on the right of figure 4.5) has different behavior but similar consequences. R$_{CHN}$ is not able to describe the whole reaction either. In the R$\to$S direction, after His297 abstracts the hydrogen from the substrate (R $_{CHN}\sim-1$) and after some progress in the configuration inversion, there is a simulation window where the molecular dynamics falls into the steep valley directly to the reactant S region without describing the other half of the reaction. This different behavior with the R$_{NHC}$ case can be attributed to the fact that the central TS (the saddle point located in chapter 3) is rather product-like (R-like). Consequently the reaction coordinate R$_{CHN}$ is able to reach some configurations that belong to the TS and falling down to the bottom at the S configuration.

The two flat zones encountered after the two corresponding proton transfers (in figure 4.5 are around R $_{NHC}\sim$ 0.4 and R $_{CHN}\sim$ -0.5) and before the configuration inversions step are not thermically stable at 298 K. Therefore, although these two structures were minima with very low barrier saddle points, they cannot be considered intermediates in terms of free energy. As a consequence, the mechanism cannot be considered stepwise and a unique reaction coordinate should be able to describe the whole step.

The energetics for the free energy profiles are collected in the table 4.3. Despite the highest energy point is meaningless, the energy corresponding to the inflection zones is well described and they will be valuable in the following sections.

Table 4.3: Energetics in kcal/mol corresponding to the PMF using R$_{NHC}$ and R$_{CHN}$. In brackets the corresponding reaction coordinate bin ($\pm$0.005 Å)

Reaction coordinate	R$_{NHC}$	R$_{CHN}$
Minimum	S: 0 (-0.675)	R: 0 (0.575)
Inflection barrier	14.29 (0.305)	15.59 (-0.405)
Inflection minimum	14.25 (0.365)	15.53 (-0.425)
Highest point	19.46 (1.075)	19.44 (-1.045)

We tried to include the coordinate r$_\alpha$ (see table 4.1) to describe the central step where the configuration inversion takes place. However, any movement of the substrate substituents bounded to the alpha carbon (hydroxyl, phenyl and carboxyl groups) provokes a remarkable variation in the value of $r_\alpha$ and this makes the scanning of this coordinate a very difficult task.

HCH:
Another option for Rc still using a combination of two distances consists in employing the distance between the two transferred hydrogens with respect to the alpha carbon R$_{HCH}$=r$_{HC}$-r$_{CH}$.

In figure 4.6 we show the PMF and the P(Rc) for every simulation window using R$_{HCH}$. In some regions the restraining harmonic constant was increased to obtain a better sampling and an adequate overlap between consecutive windows. The PMF profile shown in the left of figure 4.6 has been built by WHAM technique. The PMF profile by matching the different windows (not shown) had some sharp points in the regions where WHAM displays a discontinuity.

The discontinuity in the free energy profile shows clearly the inadequacy of this reaction coordinate. Although we have been able to connect the reactants and products through different overlapping simulations, this reaction coordinate describes the racemization process but lacks of the adequate exploration of some regions along the reaction process.

Figure 4.6: Left: The PMF profile using R$_{HCH}$ as reaction coordinate. Right: The histogram showing P(Rc) for the different simulations. In this particular case some windows had a bigger restriction umbrella potential to improve the sampling

This fact means that there are relevant movements that belong to the reaction path not contained in R$_{HCH}$ when scanning the distances H $_{166}\cdots$C and C$\cdots$H$_{297}$. A confirmation of this assessment is shown in figure 4.7 where the average and the fluctuations of the two bond distances r$_{HC}$ and r$_{CH}$ for each simulation window are shown. The evolution of the distances shows an abrupt jump that explains the discontinuity in the free energy profile. In particular, from S$\to$R direction, the initial approximation of Lys166 to the substrate needed for the first proton transfer is not reproduced and the proton transfer takes place at Rc $\sim$ -0.8 too brusquely. This sudden change is represented in the PMF profile by the discontinuity. The situation is the same in the R-side for the proton transfer between His297 and the substrate. Although the configuration inversion is monitored adequately, the posterior proton transfer from His297 takes place again with a brusque jump between two adjacent and overlapping windows around Rc $\sim$ 0.6.

Figure 4.7: Evolution of the H $_{166}\cdots$C (left) and C$\cdots$H$_{297}$ (right) distances. The average value for every window is plotted along with the standard deviation.