Combining four bond distances

Since the precedent two-distances reaction coordinates failed in the PMF computation, we included the four bond distances that describe the two proton transfers, namely:

$$R_4\equiv R_{NHC-CHN}=R_{NHC}+ R_{CHN}=r_{HC}-r_{NH} + r_{HN}-r_{CH}$$ (5.2)

Intuitively, we thought that if R$_{NHC}$ reproduced adequately the S-side of the reaction and R$_{CHN}$ the R-side, it indicates that joining the two of them could be an adequate way to describe the whole process.

In addition, this coordinate gives freedom to the variation of the four distances since the biasing potential only penalizes the whole combination R$_4$. Therefore this would permit the observation of an hypothetical stepwise mechanism. R$_4$ is scanned from negative values (see table 4.2) by increments of 0.2 Å until $\sim$2.5 Å.

In figure 4.8 the histogram of P(Rc) and the PMF profiles are shown. The two free energy profile correspond to the one obtained by WHAM technique and the other by matching the different windows.

Figure 4.8: Left: the PMF profile using R$_4$ reaction coordinate. The two profiles correspond to the PMF built by WHAM or by a direct match of the overlapping windows.Right: Histogram diagram displaying the probability distribution.

Despite of some difficulties that we comment below referring to the differences in the statistical analysis of P(Rc), the PMF computed using this reaction coordinate is able to draw the free energy profile all along the process. The shape of the profile is not the same as the one we encountered using R$_{NHC}$ in the S-side and R$_{CHN}$ in the R-side (figure 4.5), the main difference is that the inflection zone in the R-side is now less detailed. However, the energy for the early small barriers representing the proton transfer is approximately the same.

In addition, and what makes R$_{4}$ the most adequate coordinate, is that now we can describe the central step, where the configuration inversion takes places and which is the highest point in the PMF, along with the approximation of Lys166 and His297 for the two proton transfers.

PMF: WHAM vs matching:
There are some differences between the results coming from WHAM analysis and by adjusting automatically the several adjacent windows by a simple match criterion. The free energy profile is in some points up to $\sim$1 kcal/mol different. As we already said in section 1.5.2, WHAM technique is a more sophisticated iterative procedure that takes into account all the data without discarding the overlapping regions in the simulation and that avoids the uncertainty involved in the matching process. However, we cannot ensure for certain that WHAM method gives the correct energy without a deeper analysis. In order to examine which profile is the right one, a comparison between the two free energy profiles obtained for R$_4$ with those obtained previously will be useful. In table 4.4 the energetics of the different points is shown. Since R$_4$ is unable to describe the R-side as accurately as R$_{CHN}$ we only give an inflection point as the point where the two slopes cross. A comparison between table 4.4 and table 4.3 can be made at the two inflection points. If we accept that R$_{NHC}$ and R$_{CHN}$ can reproduce adequately the S-side and R-side respectively we must conclude that the better free energy profile computed with R$_4$ will come from the matching technique. The energy corresponding to the inflection points in table 4.3 are very similar to the points Sts, Is and Rts in table 4.4 matching the windows.

Table 4.4: Free energy in kcal/mol for the PMF using R$_4$ as a reaction coordinate. In brackets the corresponding reaction coordinate bin ($\pm$0.005 Å)

PMF Technique	Matching	WHAM
Minimum S	0 (-3.095)	0 (-0.395)
Inflection Barrier (Sts)	14.202 (-1.395)	15.25 (-1.475)
Inflection Minimum (Is)	14.116 (-1.355)	15.20 (-1.415)
TS	19.38 (0.575)	20.70 (0.455)
Inflection (Rts)	18.13 (1.135)	19.00 (1.125)
R	2.32 (2.675)	1.49 (2.675)

In figure 4.9 a detail of the matching windows process is shown. The different free energy profiles are matched at the maximum overlap point. What is shown in figure 4.9 is the unbiased and matched free energy profiles for the different windows. In this case it is included the overlap zones that will be discarded to obtain a single free energy profile (left of figure 4.8). We can see that there exist a simulation at the left of the graphic 4.9 in which the free energy increases abnormally with respect to the others. This simulation window is exploring a bad region in the R$\to$S direction before the hydrogen from His297 is transferred. This fact explains the difference between the two PMF profiles computed by WHAM and by direct matching of the windows. While WHAM takes into account all the given data, including this abnormal simulation, the matching process only employs the most populated bins and discard the overlapping regions. In a posterior examination we discarded the abnormal simulations in the WHAM analysis and then the free energy profiles from the two techniques coincided and both gave the same free energy profile.

Figure 4.9: A detail of the matching windows process. Free energy profiles for every window without removing the overlapping region.

In the histogram of figure 4.8 we can see how the probability distribution of the reaction coordinate between two consecutive windows have a good overlap. However, the difference between the WHAM analysis and direct matching of windows forces us to think that there must be an abrupt change in some of the distances involved in the reaction coordinate.

A further analysis comes from monitoring the evolution of the four bond distances during the reaction. The plots displayed in figure 4.10 show the average of the four distances that combine in R$_4$ at every simulation window.

Figure 4.10: Evolution of the four bond distances in R$_4$. The average value for every window is plotted along with its standard deviation.

The distances r$_{HC}$ and r$_{CH}$ do not change brusquely, but the approximation of Lys166 and His297 represented by the distances r$_{NH}$ and r$_{HN}$, respectively, is in some zones not well reproduced. An analysis of the geometries obtained during the reaction shows that when the approximation of one residue is not energetically favorable it is compensated in the reaction coordinate combination by the approximation of the conjugated residue whose movement is more labile. For example, when the first proton abstraction by Lys166 to the substrate must take place, His297 which is not coordinated to any residue and whose movement is rather free, tends to approximate to the (S)-substrate in order to compensate in the four distances combination R$_{4}$ the more energetic proton abstraction by Lys166. After the proton transfer has taken place His297 goes back to its position.

This situation is repeated more remarkably in the R-side with the Lys166 free movement in order to compensate the proton transfer between the (R)-mandelate substrate and His297. This is why residue Lys166 and His297 go forwards and backwards during the corresponding proton transfer in which they do not participate actively.

A solution to this problem could be weighting the distance combination in order to give preference to the important movements. For example, during the S-side of the reaction where Lys166 must abstract the proton, we could penalize the movements of His297 and do the same for the R-side with Lys166. We tried different weighting schemes mainly based on the evolution of the four distances obtained by optimization techniques (table 4.1) but the results obtained did not improve the sampling already obtained without the weighting.