Hydrogen bond analysis between a protein and a small molecule
You want to study the hydrogen bonds between a small molecule and a protein.
We can use the hbond command that will use a geometric consideration (distance of 3.0 Å and an angle cutoff of 135°) to consider a hydrogen bond is formed. For this recipe we are going to use a 1 microsecond simulation of the human papilloma virus E6 protein and the natural product kappa-carrageenan. A short trajectory of the system is presented in the following video:
We can see from the trajectory that the ligand (using spacefill representation) is interacting with the binding site of the protein throughout the simulation. We will use the topology file vph.prmtop and the trajectory file vph.nc. The system consists of a protein receptor with 151 amino acids and two zinc ions. The ligand residue id is 154. The trajectory consists of 5074 frames and has been pre-processed to remove the waters and re-oriented with the autoimage command.
parm vph.prmtop trajin vph.nc 1 last 1 hbond contacts :1-154 avgout avg.dat series uuseries hbond.gnu nointramol go lifetime contacts[solutehb] out contacts-lifetime.dat go
The CPPTRAJ scripts reads the topology file and the trajectory file using every frame. Then we call the hbond command that will have the name ‘contacts’. The command will look for Hbond donors/acceptors in the regions specified by the mask that includes residues one through 154. Very important here is the nointramol keyword, which will make the hbond command to only look for hbonds in different molecules. If we forget the keyword, the command will look for hbonds between all the 154 residues. The avgout keyword will write solute-solute hydrogen bond averages to a file called avg.dat. The series keyword will save hydrogen bond formed (1) or not formed (0) per frame for any detected hydrogen bond. The uuseries keyword will write the raw hydrogen bond time series data to hbond.gnu, which can be visualized using gnuplot. The nointramol keyword will make sure that only hydrogen bonds between different molecules are calculated.
We have to run the command to actually look for hydrogen bonds between the protein and the ligand using the go keyword. As an extra analysis, we will calculate the lifetime of those hydrogen bonds found. The lifetime will use as input the previously calculated data from the hbond command and print the results to contacts-lifetime.dat
Running the CPPTRAJ script will generate the files: avg.dat, hbond.gnu, contacts-lifetime.dat and crv.contacts-lifetime.dat. The avg.dat file will print the average of each solute-solute hydrogen bond (sorted by population) formed over the course of the trajectory. Looking at the first 10 lines of the avg.dat, we have:
#Acceptor DonorH Donor Frames Frac AvgDist AvgAng CYS_51@O LIG_154@HO2 LIG_154@O20 4536 0.8940 2.7367 163.5635 LIG_154@O19 TYR_32@HH TYR_32@OH 3686 0.7264 2.7974 163.7224 LIG_154@O22 CYS_51@H CYS_51@N 1367 0.2694 2.8457 162.6448 LIG_154@O24 CYS_51@H CYS_51@N 1243 0.2450 2.8504 162.6914 LIG_154@O7 TYR_54@H TYR_54@N 882 0.1738 2.9019 161.8689 LIG_154@O23 CYS_51@H CYS_51@N 872 0.1719 2.8483 162.5888 TYR_54@O LIG_154@HO LIG_154@O1 370 0.0729 2.7324 162.1249 LIG_154@O10 ARG_55@HH11 ARG_55@NH1 369 0.0727 2.8210 155.3269 LIG_154@O14 ARG_129@HH22 ARG_129@NH2 351 0.0692 2.8229 153.6910
We can see from the fraction column that for 89% of the simulation time there is a hydrogen bond between the oxygen atom of the residue CYS51 and the hydrogen atom of atom #2 of the ligand. Following this there is a hydrogen bond present for 72% of the simulation between the O19 of the ligand and the OH group of TYR32. In addition, 26% and 24% of the time of the simulation, a hydrogen bond is populated between oxygen atoms 22 and 24 of the ligand and the same CYS51 residue. The rest of the hydrogen bonds corresponds to 17% of the simulation between the ligand’s O7 and O23 oxygen atoms and the protein residues TYR54 and CYS51. From the file we can see the the rest of the hydrogen bonds are very short lived and we can assume that do not provide significant stability of the ligand-protein interaction.
The lifetime analysis allows the study of the time-dependent behavior of the hydrogen bonds. The first 10 lines of the file contacts-lifetime.dat are:
#Set lifetime lifetime[max] lifetime[avg] lifetime[frames] lifetime[name] 1 5 3 1.4000 7 LIG_154@O4-GLN_6@NE2-HE21 2 1 1 1.0000 1 LIG_154@HC7-GLN_6@NE2-HE21 3 1 1 1.0000 1 LIG_154@HC9-GLN_6@NE2-HE21 4 3 1 1.0000 3 LIG_154@O4-GLN_6@NE2-HE22 5 4 1 1.0000 4 LIG_154@HC6-GLN_6@NE2-HE22 6 1 1 1.0000 1 LIG_154@HC7-GLN_6@NE2-HE22 7 36 2 1.0556 38 LIG_154@O4-ARG_8@NH2-HH21 8 4 1 1.0000 4 LIG_154@HC7-ARG_8@NH2-HH21 9 46 2 1.0217 47 LIG_154@O4-ARG_8@NH2-HH22
Using the bash command: sort -g contacts-lifetime.dat -k5 to sort the file using the 5th column (lifetime frames) returns:
38 266 17 3.2782 872 LIG_154@O23-CYS_51@N-H
41 515 10 1.7126 882 LIG_154@O7-TYR_54@N-H
39 329 24 3.7781 1243 LIG_154@O24-CYS_51@N-H
37 383 24 3.5692 1367 LIG_154@O22-CYS_51@N-H
25 760 35 4.8500 3686 LIG_154@O19-TYR_32@OH-HH
244 276 105 16.4348 4536 CYS_51@O-LIG_154@O20-HO2
This information show that the longest consecutive number of frames containing the CYS_51@O-LIG_154@O20-HO2 hydrogen bond is 105 (in the lifetime[max] column, which still appears to be the longest). So the longest “lifetime” for this hydrogen bond is 105. It’s a subtle but important difference, since it could occur to have a hydrogen bond that is formed for a long time but maybe only 1 or 2 times (low frequency), or a hydrogen bond that is formed for most of the simulation but maybe for only 1 or 2 frames at a time (high frequency).
Another example for the Hbond command is available in this recipe