DFT calculation of the Catalytic Triad in Serine Proteases

In order to understand a charge relay we can use several types of analysis.

1. Potential energy surface for hydrogen bonds.

2. Mulliken or Hirshfeld charge on the active serine (methanol) hydroxyl group.

Calculating a Potential Energy Surface

To setup up a potential energy surface calculation you will use the hbvec (hydrogen bond vector) program for hydrogen bonds the interaction energy of the hydrogen atom with the H-bond donor and acceptor. The hbvec program specifies a vector between the donor and acceptor and moves one or more atoms along that vector to generate a series of structures. A single-point energy calculation is carried out for each of the structures and the result is plotted to generate a potential energy surface. Comparison of different environments permits an understanding of the specific effects of hydrogen bond relay and other group interactions on the activation of molecules at the enzyme active site.

For hbvec good set of displacements is 0.8 Å for the starting position and displacements of 0.05 Å. A good number is 30.

Once you have run hbvec, copy an input file and name it root1.input where root is the four-letter root name for your hydrogen bonding files. A script is also generated by hbvec that automatically copies this input file to input files for all of the other car files. The script has the name root.job. It is advisable to specify the Thermal option under the Occupation keyword. This will give a much higher probability that the self-consistent field (SCF) calculation will converge for all of the structures you have generated in hbvec.

You will also need to set the charge for this collection of molecules. The charge is -1 since the acetic acid is deprotonated in the coordinates file given to you. Using the vi editor you can search for the word Charge in the input file by typing /Charge in vi. Using the i command to enter text you can add a charge of -1. You can delete the charge of 0 (the default) using the x (delete) command.

Activate the script generated by hbvec or fragvec using the following UNIX syntax:

Ø chmod +x [script_name]

Then run the script and it will automatically copy the input file giving the appropriate names for the various displaced geometries.

Once the job is completed you can extract the energies using pesgen_thermal. Type the word ‘none’ when prompted for the second set of outmol files (since you only have one). The file you will obtain in this case is the root.energy file. The first column is the displacement. The second column is the uncorrected energy (in kJ/mol). The third column is a correction term for finite temperature. The fourth column is the corrected data (in kJ/mol). A plot of the model for the serine protease hydrogen bond number 2 is shown in Figure 1 below.

Figure 1. Potential energy surface for displacement of the hydrogen along the bond vector from the methanol (serine) oxygen to the imidazole (histidine) nitrogen.

The equilibrium bond length is almost exactly 1.0 Å. Note that the energy we have extracted is the binding energy and not the total energy. In some instances you will need to extract the total energy from the outmol file. However, most of the time we are interested in the relative energy since it is the relative energy as a function of displacement that we will want to compare for different structures. The relative energy in this case can be obtained by adding 9928 kJ/mol to all of the points in Figure 1 and is plotted in Figure 2.

Figure 2. Relative energy of displacement of the hydrogen atom along the hydrogen bond vector number 2 in the catalytic triad model.

Further comparisons can be made as indicated in the lecture and proposed in suggested research topics.

Calculating the Charge Distribution

Change the commenting on the Mulliken_Charge keyword in the input file. TO do this you can use the vi editor. Search for Mulliken_Charge using the command "/Mul". If the commenting appears as below both the Mulliken and Hirshfeld charges will appear in the outmol file when the job is completed.

# ----- Properties Keywords

#Bond_Order on

Mulliken_Analysis charge

#Mulliken_Analysis population

#Mulliken_Analysis full

Hirshfeld_Analysis charge

#Hirshfeld_Analysis dipole

#Hirshfeld_Analysis quadrupole

In order to illustrate the kind of analysis that can be performed using the Mulliken charge we will extract the charge on the methanol (serine) oxygen and on the imidazole (histidine) nitrogen as the hydrogen atom is displaced as above. Instead of using a program to do this we can use UNIX commands. Here the grep command is useful. First, you need to examine (any) one of the outmol files. Use the vi editor and search for the Mulliken charge using the “/Mul” command (backslash means search for the characters anywhere in the file, NOTE: case is important!). The first occurrence will be in the top of the file where the input file is repeated. Just type “n” to find the next occurrence. Keep typing “n” until you find a list of Mulliken charges. It will look like this.

Mulliken atomic charges:

charge spin

C( 1) -0.091 0.000

H( 2) 0.003 0.000

H( 3) 0.003 0.000

H( 4) 0.000 0.000

O( 5) -0.636 0.000

C( 6) 0.517 0.000

O( 7) -0.566 0.000

N( 8) -0.297 0.000

C( 9) 0.081 0.000

C( 10) 0.082 0.000

N( 11) -0.423 0.000

C( 12) 0.247 0.000

H( 13) 0.271 0.000

H( 14) -0.012 0.000

H( 15) -0.029 0.000

H( 16) -0.005 0.000

H( 17) 0.311 0.000

O( 18) -0.630 0.000

C( 19) 0.254 0.000

H( 20) -0.025 0.000

H( 21) -0.025 0.000

H( 22) -0.030 0.000

Determine which atom is the serine oxygen and which is the histidine nitrogen. In my file the oxygen is atom 18 and the nitrogen is atom 11. Now, you want to obtain a list of all of the oxygen charges in all of the files. You do not need to open 40 files to do this. Instead, exit vi using “:q” and use the UNIX “grep” command. The “grep” command finds the line that contains the text indicated. Try the following:

Ø grep ‘C( 18)’ *.outmol

A list of the line in each outmol file will flash across the screen. To capture this list to a file use the UNIX redirect as follows:

Ø grep ‘C( 18)’ *.outmol > c.charge

Now all of the charges are in the file c.charge. You can do the same for the nitrogen.

Ø grep ‘N( 11)’ *.outmol > n.charge

The lists will need to be sorted. This can be done best in vi editor.

Ø vi o.charge

The way to sort a list like this one conveniently is to cut and paste. Cut the first ten lines using the command “10dd”. The lines will disappear. Then type “p” for paste and they will reappear below the first line. Then type “20dd” and then “p” and then “30dd” and then “p”. Now the list is sorted. The list will look like this for O.

ct__1.outmol: O( 18) -0.397 0.000

ct__2.outmol: O( 18) -0.450 0.000

ct__3.outmol: O( 18) -0.498 0.000

ct__4.outmol: O( 18) -0.539 0.000

ct__5.outmol: O( 18) -0.575 0.000

ct__6.outmol: O( 18) -0.605 0.000

ct__7.outmol: O( 18) -0.630 0.000

ct__8.outmol: O( 18) -0.651 0.000

ct__9.outmol: O( 18) -0.669 0.000

ct__10.outmol: O( 18) -0.684 0.000

ct__11.outmol: O( 18) -0.697 0.000

ct__12.outmol: O( 18) -0.707 0.000

ct__13.outmol: O( 18) -0.717 0.000

ct__14.outmol: O( 18) -0.724 0.000

ct__15.outmol: O( 18) -0.731 0.000

ct__16.outmol: O( 18) -0.736 0.000

ct__17.outmol: O( 18) -0.740 0.000

ct__18.outmol: O( 18) -0.742 0.000

ct__19.outmol: O( 18) -0.743 0.000

ct__20.outmol: O( 18) -0.742 0.000

ct__21.outmol: O( 18) -0.740 0.000

ct__22.outmol: O( 18) -0.738 0.000

ct__23.outmol: O( 18) -0.735 0.000

ct__24.outmol: O( 18) -0.732 0.000

ct__25.outmol: O( 18) -0.729 0.000

ct__26.outmol: O( 18) -0.726 0.000

ct__27.outmol: O( 18) -0.724 0.000

ct__28.outmol: O( 18) -0.722 0.000

ct__29.outmol: O( 18) -0.720 0.000

ct__30.outmol: O( 18) -0.719 0.000

To manipulate this further you can read the file into Excel. When plotted as a function of distance the Mulliken charges appear as shown in Figure 3 below.

Figure 3. Mulliken charge for the three atoms involved in hydrogen bond number 2 in the catalytic triad model.

Note that as the hydrogen atom is extracted from the oxygen the Mulliken charge becomes more negative, increasing from approximately -0.6 at the equilibrium distance nearly -0.8. The nitrogen on the other hand becomes more positive as the hydrogen moves closer to it. This type of analysis can be used to rationalize trends when comparing different hydrogen bonding configurations. The simple assumption is that increased negative charge corresponds to increased nucleophilicity of the serine hydroxyl group.