Discrepancy Between Multi-Ligand Docking Score and Summed Individual Scores in AutoVina

[Squirrelity]

I am using AutoVina for multi-ligand docking and encountered an issue related to the binding energy scores. A fter performing a multi-ligand docking run, I extracted each ligand from the complex and calculated their binding scores individually using the "score only" mode. I observed that the sum of these individual scores is consistently higher than the binding energy score of the original multi-ligand docked complex.

I calculated multi-ligand docking score by these code:
`
rep=rep
lig=[A,B]

v = Vina()
v.set_receptor(rigid_pdbqt_filename=rep)
v.set_ligand_from_file(lig)
v.compute_vina_maps(center=center, box_size=box_size)
v.dock(min_rmsd=0.5)
v.dock(exhaustiveness=40)
v.write_poses(pdbqt_filename=out_dir+"%s_dockingPose.pdbqt"%("_".join(lig_name)))
print("file : %s_dockingPose.pdbqt is saved."%("and".join(lig_name)))
`

I calculated ligand docking score by these code:
`

docking ligand1

v = Vina(no_refine=True)
v.set_receptor(rep)
v.set_ligand_from_file(lig[0])
v.compute_vina_maps(center=center, box_size=box_size)
score = v.score()
print(f"Binding affinity:{lig_name} {sum(score)} kcal/mol")

docking ligand2

v = Vina(no_refine=True)
v.set_receptor(rep)
v.set_ligand_from_file(lig[1])
v.compute_vina_maps(center=center, box_size=box_size)

score = v.score()

print(f"Binding affinity:{lig_name}_SAM: {sum(score)} kcal/mol")
`

I use pymol to get the center and box size of the ligand.

Result:
The multi-ligand docking score:(A+B)
REMARK VINA RESULT: -11.001 0.000 0.000
and each ligand docking score is follows:
Binding affinity:A: -19.916 kcal/mol
Binding affinity:B: -8.59 kcal/mol

Objective:
I aim to understand the contribution of each ligand to the overall binding affinity in a multi-ligand context. However, the summed individual binding scores do not match the binding score of the multi-ligand complex, leading to confusion about how to interpret the results.

Questions:

1. Proper Methodology: What is the recommended approach to calculate and interpret the contribution of each ligand in a multi-ligand docking scenario using AutoVina?
2. Scoring and Energy Calculation: Are there specific settings or methodologies in AutoVina that should be used to achieve accurate individual ligand contributions from multi-ligand docking results?
3. Common Pitfalls : Are there known issues or adjustments to be aware of when comparing individual ligand scores to the multi-ligand docking score?

AutoVina version:

[rwxayheee]

Hi @Squirrelity
This is because of how the number of torsions is counted for a multi-fragment ligand. I encourage you to check the following:

In your dual-ligand output, you can retrieve the number of torsions by the following values printed after REMARK:
total torsion numbers = (INTER / VINA RESULT - 1)/0.0585

From singly scored poses, you can retrieve the number of torsions by the following items in v.score():
individual torsion numbers = (score[1]/ score[0] - 1)/0.0585

The individual torsion numbers should add up to the total torsion number. As you can see, the multi-ligand poses are penalized more (scaled down) because of the degrees of freedom.

I'm not very sure about the best scoring function for docking multiple ligands. But one common pitfall i can think of is that the Vina scoring function penalizes ligands with higher degrees of freedom by a scaling factor. Therefore, I would avoid comparing ligands that are drastically different in terms of size and/or number of free torsions. The scoring function might be more useful (and fair) when comparing ligands or multi-ligand systems that are similar in size and number of free torsions.

Please let us know what you think and if you find any further issues with that!

[Squirrelity]

Hi @rwxayheee
Thank you for the detailed explanation regarding the discrepancy due to calculations in multi-ligand docking.

I have followed your instauctions to calculate the individual torsion numbers and the total torsion number. Here are my results:

Individual Torsion Numbers:

For ligand A: 6.996327125672311 kcal/mol
For ligand B: 6.996979628558575 kcal/mol

Calculation Details for Total Torsion Numbers from Multi-Ligand Output:
MODEL 1 REMARK VINA RESULT: -11.001 0.000 0.000 REMARK INTER + INTRA: -22.274 REMARK INTER: -20.004 REMARK INTRA: -2.270 REMARK UNBOUND: -2.270

Using the formula (INTER / VINA RESULT - 1) / 0.0585 ,I computed the total torsion number to be 28.49335204972235, which does not match the sum of the individual torsion numbers.

This discrepancy between the total torsion number and the sum of individual torsion numbers is causing confusion. Could you please clarify whiy the total number from the multi-ligand output does not math the sum of the individual torsion number?

[rwxayheee]

Hi @Squirrelity
Thanks for the information. Could you please share an example of your docking commands, and receptor & ligands files?

I thought the number of torsions added up after checking my own example of docking two copies of the same ligand, but maybe I'm wrong.. thanks for reporting, and I can take a closer look at your files if that helps

[Squirrelity]

Hi @rwxayheee
My ligand is v22 from 7O2I（pdb）and sam from 5K7U(pdb)，my receptor is 5K7U.
I used python api to dock. My docking commands are as follows:
`lig_list = [[v22],[v22,sam]]
for lig in lig_list:
rep = 5K7U
lig_name = [i.split("/")[-1].split(".")[0] for i in lig]

v = Vina()
v.set_receptor(rigid_pdbqt_filename=rep)
v.set_ligand_from_file(lig)
v.compute_vina_maps(center=center, box_size=box_size)
v.dock(min_rmsd=0.5)
v.dock(exhaustiveness=40)
v.write_poses(pdbqt_filename=out_dir+"%s_dockingPose.pdbqt"%("_".join(lig_name)))
print("file : %s_dockingPose.pdbqt is saved."%("_and_".join(lig_name)))`

I found the "INTER" and "Vina result" from dockingPose.pdbqt

and I calculate individual torsion number by using this code:
v = Vina(no_refine=True) v.set_receptor(METTL3_METTL14_5K7U) v.set_ligand_from_file(lig) v.compute_vina_maps(center=center, box_size=box_size) score = v.score() print(f"Binding affinity:{lig_name}_SAM: {sum(score)} kcal/mol") individual_torsion_numbers = (score[1] / score[0] - 1)/0.0585 print(f"individual_torsion_numbers:{individual_torsion_numbers} kcal/mol")

[rwxayheee]

Hi @Squirrelity
I will set up a similar calculation in a moment. After taking a quick look at the PDB structures, one thing I noticed in your original post is that -
When you were doing "docking ligand1" and "docking ligand1", you didn't v.dock(). Was that intentional? If you didn't ask Vina to dock, the scores will be from single-point evaluation. The original positions of the ligands in their inputs would not change, and the scores likely won't add up to a co-docking score, especially if the original positions overlap
I will write another reply with my docking results shortly

[rwxayheee]

Hi @Squirrelity

With just v.dock(), the torsions should add up in these single point calculations. Here's what I did & what I got:

5K7U_A.pdbqt.txt
V22.pdbqt.txt
SAM.pdbqt.txt

from vina import Vina

v = Vina(sf_name='vina')

v.set_receptor('5K7U_A.pdbqt')

lig_V22 = 'V22.pdbqt'
lig_SAM = 'SAM.pdbqt'

bcenter = [-43.014, 34.412, 25.684]
bsize = [40, 40, 40]

v.compute_vina_maps(center=bcenter, box_size=bsize)

v.score() for V22; Ntor = 7

Input

v.set_ligand_from_file(lig_V22)

score = v.score()
print(score)
print((score[1]/score[0]-1)/0.0585)

Output

[ 32.668  46.037   0.      0.      0.     -0.651 -13.369  -0.651]
6.995528178337043

v.score() for SAM; Ntor = 9

Input

v.set_ligand_from_file(lig_SAM)

score = v.score()
print(score)
print((score[1]/score[0]-1)/0.0585)

Output

[-4.807 -7.337  0.     0.     0.    -0.634  2.529 -0.634]
8.996851102114258

v.score() for V22 and SAM; Ntor = 16

Input

v.set_ligand_from_file([lig_V22, lig_SAM])

score = v.score()
print(score)
print(((score[1]+score[2])/score[0]-1)/0.0585)

Output

[ 220.713   38.7    388.458    0.       0.      -1.285 -206.446   -1.285]
15.988973730475136

Here Vina writes the interactions between ligands in flex_inter (which is score[2]). I didn't know about this... So to calculate the Number of torsions (Ntor), we should do:
((score[1]+score[2])/score[0]-1)/0.0585?
I'm not entirely sure when other_inter (score[3]) could get a nonzero value. But in the output PDBQT, score[1] and score[2] are usually added up to REMARK INTER. Continue reading for docking calculation:

Docking V22 and SAM; Ntor = 16

Input

v = Vina(sf_name='vina',seed=666)

v.set_receptor('5K7U_A.pdbqt')
v.set_ligand_from_file([lig_V22, lig_SAM])
v.compute_vina_maps(center=bcenter, box_size=bsize)

v.dock()
score = v.score()
print(score)
print(((score[1]+score[2])/score[0]-1)/0.0585)

Output

[-11.422 -21.994  -0.111   0.      0.     -1.607  10.683  -1.607]
15.988039276430099

PDBQT Output

MODEL 1
REMARK VINA RESULT:   -11.422      0.000      0.000
REMARK INTER + INTRA:         -23.712
REMARK INTER:                 -22.105
REMARK INTRA:                  -1.607
REMARK UNBOUND:                -1.607

((-22.105/-11.422)-1)/0.0585 = 15.9880392764 ... approx. 16 = 7+9
Again,VINA RESULT is not a direct summation of individual Vina scores because of the scaling by number of torsions. However, the additive relation of individual torsions (16 = 7+9) still stands to me. But please so let us know if you find any further issues with anything. Thanks a lot for bringing this up