Molecular docking

AutoDock Vina predicts how a small molecule, the ligand, binds to a protein, the receptor. It is interfaced by the VinaApp in the biotite.application_v2.autodock subpackage. Let’s dock biotin into its famous binding partner streptavidin. As we would like to check how well the docking works, we download a high-resolution crystal structure of the complex, which we split into the receptor and the reference binding pose of the ligand. The docking requires bonds and formal charges, hence both are included when loading the structure.

from tempfile import gettempdir
import biotite.database.rcsb as rcsb
import biotite.structure as struc
import biotite.structure.io.pdbx as pdbx

pdbx_file = pdbx.BinaryCIFFile.read(
    rcsb.fetch("2rtg", "bcif", gettempdir())
)
structure = pdbx.get_structure(
    pdbx_file, model=1, include_bonds=True, extra_fields=["charge"]
)
# The structure is a homodimer, one monomer is sufficient
structure = structure[structure.chain_id == "B"]
receptor = structure[struc.filter_amino_acids(structure)]
ref_ligand = structure[structure.res_name == "BTN"]

The ligand should not be taken from the complex, as it is already in its bound conformation. Instead we take the idealized model of biotin from the Chemical Component Dictionary, that comes with the required bonds and charges.

import biotite.structure.info as info

ligand = info.residue("BTN")

Now we let Vina place the ligand into the receptor. The search space is a box, defined by its center and size. Here we cheat a little bit by centering the box at the position of the reference ligand, which we would not know for an unknown complex.

import biotite.application_v2.autodock as autodock

app = autodock.VinaApp()
result = app.run(
    ligand,
    receptor,
    center=struc.centroid(ref_ligand),
    # Edge length of the box in Å
    size=[20, 20, 20],
    # A fixed seed and a single thread make the run reproducible
    seed=42,
    cpu=1,
).result()
print(result.energies)

[-6.82 -5.88 -5.81 -5.79 -5.49 -5.49 -5.49 -5.3  -5.22]

Vina reports multiple binding modes, sorted from best to worst predicted binding energy. The coordinates of the ligand in each binding mode are given by VinaResult.ligand_coord, which we combine with the input ligand to an AtomArrayStack. Note that Vina removes nonpolar hydrogen atoms, so their coordinates are NaN. As we compare only the heavy atoms with the reference pose anyway, we simply remove all hydrogen atoms.

import numpy as np

docked_ligand = struc.from_template(ligand, result.ligand_coord)
docked_ligand = docked_ligand[..., struc.filter_heavy(docked_ligand)]
ref_ligand = ref_ligand[struc.filter_heavy(ref_ligand)]
# Ensure that both models have the same atom order
docked_ligand = docked_ligand[..., info.standardize_order(docked_ligand)]
ref_ligand = ref_ligand[info.standardize_order(ref_ligand)]
# No superimposition, as the position in the binding pocket matters
rmsd = struc.rmsd(ref_ligand, docked_ligand)
for energy, deviation in zip(result.energies, rmsd):
    print(f"{energy:6.2f} kcal/mol  {deviation:5.2f} Å")

 -6.82 kcal/mol   0.44 Å
 -5.88 kcal/mol   6.75 Å
 -5.81 kcal/mol   6.36 Å
 -5.79 kcal/mol   2.37 Å
 -5.49 kcal/mol   8.03 Å
 -5.49 kcal/mol   6.34 Å
 -5.49 kcal/mol   6.33 Å
 -5.30 kcal/mol   5.44 Å
 -5.22 kcal/mol   6.41 Å

The binding mode with the best predicted energy is also the one that is closest to the experimentally determined pose. Let’s look at it visually:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6.0, 4.0))
ax.scatter(result.energies, rmsd, marker="+", color="black")
ax.set_xlabel("Binding energy (kcal/mol)")
ax.set_ylabel("RMSD to reference (Å)")
fig.tight_layout()

Vina can also treat selected side chains of the receptor as flexible during the docking. For this purpose, a boolean mask of the flexible residues is given to the flexible parameter of VinaApp.run().