In version 3.0 of AutoGrid and AutoDock , we introduced a new kind of scoring function that is used during and at the end of the dockings. It is based on the principles of QSAR (quantitative structure-activity relationships) and was parameterized using a large number of protein-inhibitor complexes for which both their structure and inhibition constants, or K i , were known. The user is encouraged to refer to the description of how this free energy function was derived in the original literature See Morris, G. M., Goodsell, D. S., Halliday, R.S., Huey, R., Hart, W. E., Belew, R. K. and Olson, A. J. (1998), "Automated Docking Using a Lamarckian Genetic Algorithm and and Empirical Binding Free Energy Function", J. Computational Chemistry, 19: 1639-1662. .
The above diagram shows the thermodynamic cycle for the binding of an enzyme, E , and an inhibitor, I , in both the solvated phase and in vacuo . Note the solvent molecules are indicated by filled circles: they tend to be ordered around the larger molecules, but when E and I bind, several solvent molecules are liberated and become disordered. This is an entropic effect and is the basis of the hydrophobic effect. The solvent ordering around E and I , when both bound and unbound, is strongly influenced by the hydrogen bonding between these molecules. These hydrogen bonds between solvent and E, and solvent and I, contribute enthalpic stabilization, and is something we can estimate in our new free energy function.
According to Hess's law of heat summation, the change in free energy between two states will be the same, no matter what the path. So we can calculate the free energy of binding in solvent by the following equation:
D G binding,solution = D G binding,vacuo + D G solvation(EI) - D G solvation(E+I)
Since we can calculate D G binding,vacuo from our docking simulation, and can estimate the free energy change upon solvation for the separate molecules E and I , and for the complex, EI , D G solvation(EI) and D G solvation(E+I) respectively, then it is also possible to calculate the free energy change upon binding of the inhibitor to the enzyme in solution, D G binding,solv . Thus, we can estimate the inhibition constant, K i , for the inhibitor, I .
A key point to bear in mind is that most parts of the new scoring function are essentially the same as the original AutoDock scoring function used in versions prior to 3.0, except that various terms in the molecular mechanics energy function have been re-scaled by new coefficients, and new terms have been introduced. These new terms include the desolvation free energy of the ligand, and an estimate of the loss of conformational degrees of freedom of the ligand upon binding.
The coefficients were derived using linear regression analysis, and we chose the linear regression model that most closely fit the observed inhibition constant data. Thus the user should not modify these coefficients lightly.
For the curious, these coefficients are defined in "
gpf3gen.awk
" and "
dpf3gen.awk
", and their variable names and values are as follows:
# # Free energy model 140n coefficients: # FE_vdW_coeff = 0.1485 FE_estat_coeff = 0.1146 FE_hbond_coeff = 0.0656 FE_tors_coeff = 0.3113 FE_desol_coeff = 0.1711
So for example, in
See Self-consistent Lennard-Jones 12-6 parameters, before multiplication by the free energy model coefficients.
and
See Self-consistent Hydrogen bonding 12-10 parameters, before multiplication by the free energy model coefficients.
, the values used in the
AutoDock
3.0 scoring function will be as given except the van der Waals coefficients and well depth energies, e, will be scaled by 0.1485, the electrostatic energy will be scaled by 0.1146, and the hydrogen bonding terms will be scaled by 0.0656. The new terms for loss of torsional degrees of freedom upon binding and the ligand desolvation free energy will be scaled by 0.3113 and 0.1711 respectively. The torsional term is actually the number of rotatable bonds in the ligand that rotate heavy atoms, multiplied by the coefficient, 0.3113. Hydroxyl rotors, for example, are not counted. This number is actually determined by AutoTors, and is written in the ligand PDBQ file after the new keyword, "
TDOF
". This can also be set in the DPF, using the keyword "
torsdof
n
0.3113
", where
n
is the number of heavy-atom rotatable bonds.
Note that the new desolvation free energy term is only calculated for aliphatic and aromatic carbon atoms in the ligand. We found that the quality of the final empirical free energy model was not affected by the inclusion of the heteroatoms N and O in this term, so these are ignored for simplicity and speed of calculation. This does mean that the ligand input PDBQ file must distinguish between carbon atoms that are aliphatic and those that are aromatic. This is done by changing the atom names of the aromatic carbons so that the initial `C' is replaced by `A'. For example, if the ligand happened to be a peptidomimetic inhibitor which contained a Phe-sidechain, then the atom names in the PDBQ file would have to be changed from CG, CD1, CD2, CE1, CE2 and CZ, into AG, AD1, AD2, AE1, AE2 and AZ, respectively. To help the user do this automatically, AutoTors has a new option that can be invoked using the `-A' flag on the command line. This will look for all planar cyclic carbons and assume these are aromatic: it will then change their atom names automatically. The user can also override the default angle AutoTors uses to determine planarity. See the Section below on AutoTors .
This also means AutoGrid will calculate two different types of carbon grid maps, one for aliphatic carbons in the ligand (*.C.map), and one for aromatic carbons in the ligand (*.A.map). These will be input by AutoDock for the docking calculations, and the internal energy parameters must also specify the values for the aromatic and aliphatic carbon atoms.
Another point to bear in mind is that
AutoGrid
3.0 can calculate smoothed pairwise potentials: the lowest energy within a user-defined distance is stored at the current position. This has the effect of `widening' the basin of affinity. This smoothing distance is set using the keyword `
smooth
' in the grid parameter input file to
AutoGrid
. It is very important that the value of `
smooth
' is not changed from 0.5Å, since the free energy function was calibrated using this setting. If smooth is set to 0.0Å, for example, the calculated free energies will probably be too high.