AutoDock Vina Gives Different Results From Different Simulations Using Same Parameters, and Different From Experimental Results

Question

[ac.research]

Hello everyone,

When I run multiple Autodock Vina independent simulations (at exhaustiveness of 100) on the SAME protein and the SAME ligand using the exact SAME parameters but with different random seeds I get different results. Each line here is the highest binding affinity of a separate docking simulation. You can see that every time I re-dock the same protein-ligand I get large differences in "the best" binding affinity (-3 to -37):

                 mode |   affinity | dist from best mode
                      | (kcal/mol) | rmsd l.b.| rmsd u.b.
                 -----+------------+----------+----------
Simulation 1        1       -10.35          0          0      random seed: 1615206230
Simulation 2        1       -10.46          0          0      random seed: -273847999
Simulation 3        1        -10.8          0          0      random seed: 2071239951
Simulation 4        1       -11.76          0          0      random seed: 1725032250
Simulation 5        1       -1.665          0          0      random seed: -172889668
Simulation 6        1       -17.23          0          0      random seed: 1076344612
Simulation 7        1       -2.826          0          0      random seed: -16921503
Simulation 8        1       -3.042          0          0      random seed: -366290784
Simulation 9        1        -3.26          0          0      random seed: -1254578750
Simulation 10       1       -3.584          0          0      random seed: -1014984749
Simulation 11       1       -3.614          0          0      random seed: 1532520213
Simulation 12       1       -3.619          0          0      random seed: 1182472761
Simulation 13       1       -3.666          0          0      random seed: -1188328766
Simulation 14       1       -36.89          0          0      random seed: 1542674337
Simulation 15       1       -37.26          0          0      random seed: 607682916

• The search space is not on the surface, it is actually a tube within a trimer (image here), you can see that the search box volume is small (1,400 angstrom^2).

• If I perform 250 separate docking simulations each at 32 exhaustiveness I get these distribution results for 3 molecules (results here).

• f we take the mode of the distribution curves we get:

MOL6 = -14.6 Kcal/mol while the paper reported -10.97 Kcal/mol (Kd = 9e-9 M)

MOL7 = -10.8 Kcal/mol while the paper reported -10.75 Kcal/mol (Kd = 1.3e-8 M)

MOL8 = -11.3 Kcal/mol while the paper reported -6.35 Kcal/mol (Kd = 2.2e-5 M)

The paper's experimental measurements were done using the Biacore. I used the following formula to convert from the paper's Kd to dG [dG=RTlnKd]:

import math
def affinity(Kd, R=0.00198720425864083, T=298):
    dG = R*T*math.log(Kd)
    print(f'dG = {dG} Kcal/mol')

Questions:

1. How do we interpret these data? what does it mean that the computed affinity is different from the experimental affinity? is it possible to consolidate the computed affinity with the experimentally determined affinity?

2. If I want to run a virtual screen to find molecules that bind with a stronger affinity than MOL6, MOL7, or MOL8, which baseline should I use, the computed affinity or the experimental affinity?