AutoDock Vina Gives Wildly Different Results From Different Simulations Using Same Parameters

Question

[ac.research]

Hello everyone,

When I run multiple Autodock Vina independent simulations (at exhaustivness 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

Which of these results should I take?

Q1. Do I take the lowest affinity energy between 15, 50, or 100 different independent simulations?

Q2. Do I run multiple independent simulations, cluster the results, and take the value of the highest cluster?

Q3. Is there an additional evaluation step to run to get a more accurate binding affinity (before an MD simulation - which takes a very long time)?

How do I interpret this type of result?

Answer 1

[Mensur Dlakic]

If a protein is large and its whole surface is used for docking, I would not expect reproducible results, especially not after 100 dockings per run. If you know the binding site, or have a pretty good idea where it is, I suggest you narrow down your docking area.

If you want to keep doing it in unbiased fashion, the last two simulations have very low energies and their numbers are similar. I would start by inspecting those two solutions.

• I would not pick the lowest-energy solution unless it is part of a cluster.
• It should be OK to combine independent runs and cluster their solutions. Given the large difference in binding affinities, chances are solid that your best cluster will be the combination of the last two runs.
• Assuming your last two runs converged to a similar binding site, it may be a good idea to run a longer simulation (1000 dockings) that is limited to that area. That way you are likely to get a more thorough sampling which will hopefully give you a closer starting point for MD.

[ac.research]

Thank you Mensur for your detailed reply.

I have spent the past 4 months actually working in this direction, let me share with you what I have reached:

1. 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).
2. I have performed the docking as follows: 250 separate simulations each at 32 exhaustiveness (results here).
3. If 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')

Clearly, MOL7 is ok because the mode of the 250 simulations reports a similar docking affinity to the paper, but MOL6 and MOL8 are not (all simulation parameters were kept constant).

Questions:

1. How do we interpret these data? specifically for MOL6 and MOL8? what does it mean that the computed affinity is different from the experimental affinity?
2. How do we consolidate the computed affinity with the experimentally determined affinity?
3. 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?

[Mensur Dlakic]

A short answer to all your questions is that I don't know what it means. You are assuming that MOL7 is correct because the numbers agree, and the other two are not because they disagree. That may or may not be true, and that goes for both cases. I have docked molecules to discern the best orientation within a known binding site, but didn't need to worry about matching predicted Kd with experiments.

The way you posted this question as a reply to a 4 months old topic makes it unlikely that anyone other than me will see it. I suggest you start a new topic with these questions as that increases your chance that others will see it as well. As a general rule, continuing a topic after a long break is less likely to yield responses compared to starting a new topic.