Small molecule binding to T4-lysozyme L99A

This is the final part of this tutorial, concerned with analysing the MD runs and constructing and interpreting the free energy curves. The first page with background information is available here, the second one with system setup information here

Running and Analyzing the MD simulations

After running all the simulations outlined in the previous page, we obtain a tremendous amount of output data. This section shows how to analyze this so that we finally obtain an estimate of the relative binding free energy of our two ligands. A first step would be to take a look a the amber mdout files, to see if anything went terribly wrong with our physical parameters. The following tables give an overview over the output files from the six different simulations:

Transformation: Remove charges from BNZ H6 - Medium: Water

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.100-5.4241.193100299.555.962.8360.70.98420.0038bnz_prod_v0_l1.out
0.200-5.7741.200100299.805.940.3361.20.98330.0033bnz_prod_v0_l2.out
0.300-5.9681.199100300.076.23-0.7358.20.98340.0042bnz_prod_v0_l3.out
0.400-6.2051.213100299.315.991.4356.70.98390.0038bnz_prod_v0_l4.out
0.500-6.4881.250100299.505.951.9347.20.98200.0035bnz_prod_v0_l5.out
0.600-6.6811.261100299.465.90-0.7354.50.98410.0037bnz_prod_v0_l6.out
0.700-6.9771.281100299.536.230.7357.20.98440.0032bnz_prod_v0_l7.out
0.800-7.2791.299100300.386.061.0354.60.98180.0043bnz_prod_v0_l8.out
0.900-7.4451.308100300.146.074.9357.40.98350.0043bnz_prod_v0_l9.out

Transformation: Transform BNZ into PHN - Medium: Water

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.1005.1883.076100300.015.940.6354.10.98240.0040bnz_prod_v0_l1.out
0.2004.4043.033100300.295.822.3354.00.98320.0037bnz_prod_v0_l2.out
0.3003.6603.010100299.675.97-0.4359.00.98350.0036bnz_prod_v0_l3.out
0.4002.5233.112100299.806.06-0.8357.10.98170.0045bnz_prod_v0_l4.out
0.5001.3673.188100300.246.09-0.1357.60.98300.0038bnz_prod_v0_l5.out
0.6000.6023.468100299.526.04-1.2362.60.98200.0036bnz_prod_v0_l6.out
0.700-0.5983.275100299.976.13-0.3359.50.98310.0037bnz_prod_v0_l7.out
0.800-1.7573.102100300.216.154.1368.20.98410.0041bnz_prod_v0_l8.out
0.900-2.6952.989100299.846.012.7357.40.98170.0037bnz_prod_v0_l9.out

Transformation: Add charges to PHN O1,H6 - Medium: Water

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.100-29.4982.153100299.786.072.3356.80.98500.0032phn_prod_v0_l1.out
0.200-30.2132.233100299.926.050.0361.60.98340.0038phn_prod_v0_l2.out
0.300-31.0812.330100300.036.080.9361.20.98350.0040phn_prod_v0_l3.out
0.400-32.1422.486100300.345.951.3350.80.98310.0033phn_prod_v0_l4.out
0.500-32.9232.598100300.495.940.1356.10.98240.0035phn_prod_v0_l5.out
0.600-34.0522.774100300.456.020.8352.70.98330.0041phn_prod_v0_l6.out
0.700-35.7683.075100300.126.171.3357.90.98260.0042phn_prod_v0_l7.out
0.800-37.6753.415100300.215.93-1.1357.90.98260.0036phn_prod_v0_l8.out
0.900-39.8733.655100300.496.063.0364.90.98310.0041phn_prod_v0_l9.out

Transformation: Remove charges from BNZ H6 - Medium: Protein

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.100-5.3800.697100299.991.670.2135.81.01260.0011t4_bnz_prod_v0_l1.out
0.200-6.3810.726100300.041.702.1134.01.01250.0012t4_bnz_prod_v0_l2.out
0.300-5.9670.823100299.871.700.9136.01.01300.0012t4_bnz_prod_v0_l3.out
0.400-6.3680.831100299.901.670.6136.81.01270.0011t4_bnz_prod_v0_l4.out
0.500-6.5020.733100300.021.681.6133.11.01310.0010t4_bnz_prod_v0_l5.out
0.600-6.8870.712100300.191.630.1137.51.01270.0013t4_bnz_prod_v0_l6.out
0.700-6.7470.764100300.091.711.1135.71.01300.0011t4_bnz_prod_v0_l7.out
0.800-6.8850.803100300.021.670.8136.01.01280.0011t4_bnz_prod_v0_l8.out
0.900-7.1040.766100300.151.681.5137.61.01310.0013t4_bnz_prod_v0_l9.out

Transformation: Transform BNZ into PHN - Medium: Protein

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.1002.2132.198100299.991.671.9135.11.01290.0012t4_bnz_prod_v0_l1.out
0.2002.5112.078100299.971.661.1134.51.01280.0012t4_bnz_prod_v0_l2.out
0.3003.5042.000100299.871.700.8134.61.01300.0011t4_bnz_prod_v0_l3.out
0.4003.9752.314100300.051.681.3134.71.01300.0011t4_bnz_prod_v0_l4.out
0.5004.1602.997100300.161.711.5137.41.01260.0014t4_bnz_prod_v0_l5.out
0.6003.3303.413100300.291.722.0133.41.01260.0011t4_bnz_prod_v0_l6.out
0.7002.9573.247100300.141.721.1135.51.01280.0013t4_bnz_prod_v0_l7.out
0.8000.9882.943100300.161.691.8133.31.01280.0013t4_bnz_prod_v0_l8.out
0.9003.4903.524100300.031.671.2135.81.01240.0011t4_bnz_prod_v0_l9.out

Transformation: Add charges to PHN O1,H6 - Medium: Protein

LambdaDVDLrmskstepsTEMPrmsPRESSrmsDensityrmsFilename
0.100-32.8311.170100300.001.691.9132.41.01290.0011t4_phn_prod_v0_l1.out
0.200-33.1641.169100300.091.660.4134.21.01260.0013t4_phn_prod_v0_l2.out
0.300-33.3851.322100300.001.720.7133.91.01260.0012t4_phn_prod_v0_l3.out
0.400-33.8281.271100299.981.741.1135.01.01270.0012t4_phn_prod_v0_l4.out
0.500-33.8401.247100300.061.690.6136.61.01280.0012t4_phn_prod_v0_l5.out
0.600-34.1641.189100300.001.681.3134.91.01330.0011t4_phn_prod_v0_l6.out
0.700-34.4541.293100299.981.691.9137.11.01300.0013t4_phn_prod_v0_l7.out
0.800-34.6591.341100300.011.701.1137.71.01270.0015t4_phn_prod_v0_l8.out
0.900-34.3881.318100299.931.671.1133.91.01310.0011t4_phn_prod_v0_l9.out

We see temperature fluctuation in the range of 5K for the ligands in water and 2K for the complexes, about what one would expect. The pressures fluctuate wildly as always when ntp=2 is used, but densities stay relatively constant. dV/dλ rms values are in the range of several kcal/mol. In an actual study, this quick glance over the results would need to be replaced with a careful analysis of all the trajectories to see if interesting conformational changes or ligand movements happen, but it shall suffice for the purpose of this tutorial. Therefore we can proceed to analysing the resulting dV/dλ vs. λ curves.

Results

The plot below shows the dV/dλ curves assembled from 9 λ points each. We can see that the curves for transformations in water are smoother than those for the ligand change in the protein, as would be expected due to the more complex conformational landscape of a protein-ligand complex as compared to an isotropic solution. The vdW-transformation in the complex shows the biggest fluctuations of dV/dλ and the free energy curve for Step 2 in the complex is the one that looks the most 'bumpy'. Since TI-free energy curves are almost always smooth curves, this can probably be attributed to insufficient sampling and running longer simulations for the transformations in complex would most likely smooth out the free energy curves considerably.

The thermodynamic effect of removing the charge on the hydrogen atom is small and very similar when done in complex or solution. Adding the charges to the phenyl group (step 3) gives much bigger changes in energy than the other two steps, but bear in mind that the result of a single transformation is not a physically significant number, only differences between the results of two transformations are. The free energy change associated with each of the steps can be calculated by numerically integrating over the ΔG(Complex-Water) curves on the right from λ=0-1. The integration can e.g. be done using this script on a list with λ, dV/dλ and rms(dV/dλ) values like that one). The trapezoid rule (values for λ=0 and 1 were linearly extrapolated from the closest two values) gives:

Step Free Energy changerms

ΔG (Step 1)0.051.46
ΔG (Step 2)1.724.19
ΔG (Step 3)0.063.04

ΔG (Total) 1.835.34

Interestingly, the main contribution to the binding free energy difference stems from vdW interactions and not from electrostatics. One would have expected that the loss of solvation on the polar hydroxyl group might strongly influence the relative binding free energies, but from this calculations, the main difference seems to be the sterical effect of accomodating the bulkier hydroxyl group in the binding site, with electrostatics playing a minor role. This result might however change if the calculation would be refined further.
Our prediction is that phenol binds 1.8 kcal/mol weaker to T4-lysozyme than benzene. While this is not in perfect agreement with experimental data, it captures the trend of benzene being the stronger binder. To improve the result, more simulations could be added at different λ points (this is indeed one of the strong points of TI, you can add as many additional data points as you want to refine your result without having to redo the initial calculations) or more sophisticated numerical integration schemes could be used. The rms value of the final result should not be mistaken for the statistical uncertainty of the result. To get a good error estimation, some additional analysis is required (see e.g. Steinbrecher T. et al., J Chem. Phys. 127, 214108, 2008).


This tutorial was written by Thomas Steinbrecher, TSRI 2007