Chorismate mutase free Energy of reaction using UMbrella sampling

We will continue with the research we are conducting on the enzyme chorismate mutase and catalysis of the conversion of chorismate into prephenate. In this case, our objective is to calculate the free energy of the reaction for this system (different from the previous step, where we obtained a potential energy profile of the process). We will use a powerful technique known as Umbrella Sampling and for this we will need the data and coordinates obtained previously; Therefore, it is important that, if you have not already done so, you visit the tutorial on how to sample the reaction potential energy profile for this system.

The pDynamo serialization file for this system can be downloaded here.

The file containing the coordinates obtained by scanning the reaction coordinates in the previous step can be obtained here.

To access the window for calculations using umbrella sampling, simply click on the icon shown in Figure 1, on the main window toolbar, or access via the main menu > Simulation > Umbrella Sampling

The strategy we will adopt is to use each of the coordinates (frames) obtained in the previous tutorial as a starting point for a restrained molecular dynamics (MD) simulation using a harmonic potential. The restraining harmonic potential is applied according to the chosen reaction coordinate, similar to what happens for reaction coordinate scans.

The advantage of using the coordinates obtained previously is that now this computationally expensive process can be done in parallel, where each coordinate triggers an independent process. Each of these processes is called a "window."

To perform this calculation in parallel, the user should choose the 'From Trajectory' option and then specify the folder containing the set of coordinates that will serve as a reference for the calculations. The number of CPUs should be determined based on the availability of each one; in this case, 10 CPUs were used.

The reaction coordinate for this simulation is the same as the one selected for the previous tutorial, as shown in Figure 2.

In this tutorial, we will conduct a sampling of 10 picoseconds to equilibrate the system, followed by another 20 picoseconds in what we term 'data collections.' The latter period will be utilized to calculate the potential of mean force (PMF) and, consequently, the free energy of the reaction. It is important to note that we opted for a small sample size in these two stages to expedite the process. For research purposes, it is advisable to extend the sampling time. Currently, the literature recommends times of 20 picoseconds for equilibrium and 50 picoseconds for data collection.

This calculation is expected to take a few hours. As a result, a new folder will be generated containing a log file and two subfolders: one for equilibrium sampling and another for data collection. To derive the potential of mean force (PMF), it is necessary to submit the results of the data collection phase to a Weighted Histogram Analysis Method (WHAM). To do this, access the WHAM window, in main menu>analysis> WHAM. In the WHAM window, you must import all folders obtained during data collection.

The WHAM analysis results in two output files. The first, named as per user stipulation (in this case, called output.log), encompasses all the information regarding the calculated PMF. Additionally, another file is generated and appended with the 'pmf' tag, containing the PMF data—likely the information you'll need for inclusion in your scientific manuscript. Figure 5 illustrates the obtained PMF in this tutorial and how this structure is presented in the log file.

In this tutorial, we acquired a free energy profile for the reaction catalyzed by the enzyme chorismate mutase. The obtained energy barrier value is significantly higher compared to the experimental value (almost twice the expected valor). This discrepancy can be attributed to two factors: firstly, the use of the semiempirical Hamiltonian PM6, which lacks high accuracy, and secondly, the selection of a small Quantum Mechanics (QM) region to reduce sampling time. We plan to revisit this system in the future for more refined calculations aiming for truly accurate results.