Chorismate Mutase Reaction using Nudged elastic band
Introduction
In this tutorial, we will delve into the application of the Nudged Elastic Band (NEB) and QC/MM to simulate the reaction catalyzed by the enzyme chorismate mutase (CM). This reaction is characterized by a one-step conversion from chorismate to prephenate, as illustrated in Figure 1. The experimental energy barrier for this process is approximately 62 kJ/mol. It is noteworthy that this reaction does not involve any covalent intermediate with the enzyme. Instead, the enzyme primarily functions to guide and stabilize the substrate toward the near-attack conformation.
This tutorial assumes that you are proficient in pruning a system using pDynamo/EasyHybrid and have the knowledge to accurately define and set up the Quantum Region for a hybrid QC/MM system.
Figure 1: Conversion of chorismate into prephenate. Adapted from Jianpeng and collaborators (https://www.pnas.org/doi/10.1073/pnas.95.25.14640)
System Loading
First, load the system containing the CM enzyme and substrate into EasyHybrid (Figure 2). The previously prepared system is available in the file provided below:
To open this file, navigate to: Main Menu > File > Open (Change the file type from ".easy" to ".pkl"). For more details check the user guide.
Figure 2: Enzyme + substrate complex loaded into EasyHybrid. The region previously determined to be quantum is presented in the form of balls and sticks.
The system was constructed using the amber99SB force field and comprises 8061 atoms, with 5673 being frozen (atoms located outside a 14-angstrom sphere centered on the C5 carbon atom of the substrate). The quantum region consists of 24 atoms (the chorismate molecule) and is described by the semi-empirical Hamiltonian PM6. For the sake of expeditious completion of the tutorial, we opted for the smallest possible quantum region.
Harmonic Distance Restrictions / Sampling The products
The system initially containing the enzyme-substrate complex was already optimized (meaning, nothing needs to be done regarding the reactants' condition). On the other hand, we need to obtain a set of coordinates related to the products, a task that can be done in several different ways. In this tutorial, we will obtain the coordinates of the products from the coordinates of the reactants, applying a harmonic distance restraint.
This is done by selecting the C1 and C9 atoms of the chorismate molecule using the "picking selection" model and adding a harmonic restraint through the window accessed in the glmenu (Figure 3).
Figure 3: Applying the harmonic distance constraint. It is important to define an equilibrium distance close to the chemical bond you want to form.
Select a distance that is consistent with the equilibrium bond length for the formed bond. In this case, a reasonable value is 1.5. When a distance constraint is applied (and active), it is visually represented as a dashed line overlapping the atom representation. You can enable or disable the imposed restriction by accessing the selections and restrictions window.
Geometry Optimization
Proceed with 600 steps of geometry optimization using the default configuration of EasyHybrid.
Figure 4: coordinates before (left) and after (right) the geometry optimization with the harmonic constraint.
Now that we have obtained an approximate structure for the product, it is necessary to deactivate the harmonic constraint and reoptimize the geometry. This step should be done quickly, and the result will be a structure practically identical to the previous one, with just a slight deformation in the C1-C9 connection distance (data not showed).
Applying of NEB Method
With the structures of the reactants and products properly sampled, we will now apply the NEB method to simulate the reaction. It's worth noting that this method offers a significant advantage over the reaction coordinate scans discussed in other tutorials. In NEB, there's no need to define a reaction coordinate; you only need the coordinates of the reactants and products. To access the NEB window, click on the icon shown in Figure 5, located on the toolbar of the main window.
Figure 5: Icon to access the NEB tool in the EasyHybrid graphical environment. The same window can be accessed through the main menu by going to: simulation> Nudged Elastic Band (NEB).
The settings employed in this tutorial are depicted in Figure 6. If you decide to revisit this tutorial later, it is recommended to experiment with certain parameters to observe their impact on the results, such as spring force, the number of structures used, and the maximum number of iterations. For this tutorial, we will designate the final and initial images as fixed.
Figure 6: Window to run NEB calculations. Only the initial and final coordinates of the reaction are needed to execute this routine. In this case, "Reactants" are the initial coordinates, previously optimized, and "Products" are the last set of coordinates obtained.
The outcome of this simulation is a folder containing coordinate files in pkl format, accompanied by a log file named 'output.log.' This log file encompasses all pertinent information regarding the simulation. Towards the end of the file, there is a table displaying potential energy values for each set of coordinates, referred to as images. From this table, we can derive essential thermodynamic properties, including the energy barrier.
----------------------------------------------------------------------------------------------------------------------------------------------------------
Path Summary
----------------------------------------------------------------------------------------------------------------------------------------------------------
Image Abs. Function Rel. Function RMS Gradient Angle Curvature Distance (i,i+1) Distance (i,i+2) Total Distance
----------------------------------------------------------------------------------------------------------------------------------------------------------
0 -24404.519 116.860 0.000 0.0 0.000 0.586 1.149 0.000
1 -24404.741 116.638 0.050 155.5 0.000 0.589 1.165 0.586
2 -24404.277 117.102 0.026 161.6 0.000 0.591 1.145 1.175
3 -24403.205 118.173 0.047 152.4 0.000 0.588 1.122 1.766
4 -24401.208 120.171 0.042 145.8 0.000 0.586 1.131 2.354
5 -24398.614 122.764 0.077 147.8 0.000 0.591 1.136 2.940
6 -24401.376 120.003 0.079 146.5 0.000 0.595 1.119 3.531
7 -24406.525 114.853 0.055 140.7 0.000 0.593 1.149 4.126
8 -24408.821 112.558 0.071 152.9 0.000 0.589 1.011 4.719
9 -24410.078 111.301 0.061 117.9 0.000 0.591 1.100 5.308
10 -24410.502 110.876 0.092 135.7 0.000 0.596 1.117 5.899
11 -24406.586 114.793 0.157 139.0 0.000 0.597 1.043 6.495
12 -24383.201 138.178 0.285 122.6 0.000 0.592 1.061 7.092
13 -24322.840 198.539 0.241 127.9 0.000 0.589 1.134 7.684
14 -24443.451 77.928 0.367 147.4 0.000 0.593 1.150 8.273
15 -24483.463 37.915 0.262 151.7 0.000 0.593 1.132 8.865
16 -24500.738 20.640 0.126 147.4 0.000 0.587 1.101 9.459
17 -24506.630 14.749 0.131 141.9 0.000 0.578 1.119 10.045
18 -24509.820 11.559 0.045 150.8 0.000 0.578 1.151 10.624
19 -24513.423 7.956 0.048 167.3 0.000 0.580 1.123 11.201
20 -24517.614 3.765 0.036 151.3 0.000 0.579 1.138 11.781
21 -24519.880 1.499 0.038 159.7 0.000 0.577 1.153 12.360
22 -24520.849 0.530 0.026 172.9 0.000 0.579 1.142 12.938
23 -24521.379 0.000 0.046 159.4 0.000 0.582 0.000 13.516
24 -24521.314 0.065 0.000 0.0 0.000 0.000 0.000 14.098
----------------------------------------------------------------------------------------------------------------------------------------------------------
Of particular interest in this case are the relative energies (Rel. Function), indicating energy values of 116,860 kJ/mol for the reactants, 198,539 kJ/mol for the transition states, and, of course, 0.00 kJ/mol for the products. Consequently, there exists a barrier of 81.679 kJ/mol (or 19.54 kcal/mol), not bad for a simple QC/MM model. For a visual inspection of the obtained results, simply import them as suggested in Figure 7.
Figure 7: To view the results obtained, simply go to File>Import Data... Select the format "pkl folder - pDynamo trajectory", folder containing the obtained trajectory and the log file (this is optional). also define the name of the new object that will be created (here called "Reaction").
Explore the interactive graph displaying energies obtained by the NEB for each image/frame, linked with the trajectory frames, in the Energy Surface Analysis window (Figure 8).
Figure 8: Visualization of the results obtained by NEB simulation in EasyHybrid.