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.

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:

chorimate mutase pkl file.

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.

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).

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.

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.

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.

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.

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).