Flexible Vina Docking for Cyclodextrins

[jaketanderson]

Hi @diogomart,

I'm trying to optimize Vina parameters to get better agreement with experiment for docked binding energies of small molecules to cyclodextrins. To do this, I'd like to be able to take a cyclodextrin derivative and keep its main sugar ring rigid but customize the flexibility of substitutions along the ring. I'm able to use obabel to successfully dock small molecules to a rigid 6-phenyl substituted $\beta$-cyclodextrin but I need help in setting up flexibility. There is a flexible receptor setup script at example/autodock_scripts/prepare_flexreceptor.py, but I think this only works for proteins.

The host and ligand files I am using can be found in this gist. The file called bcd_deriv_flex.pdbqt was an attempt at creating a flexible 6-phenyl substitution via obabel translation of the substitution's isolated pdb to rotatable pdb. It did not work; it didn't give an error, just crashed my Jupyter notebook when I tried to run it.

The code I am using to test docking these structures is this:

from vina import Vina

v = Vina(
    sf_name='vina',
    seed=42, 
    cpu=1,
    verbosity=2
)

# weights = [-0.035579, -0.005156,  0.840245, -0.035069, -0.587439, 0.05846]
# v.set_weights(weights)

v.set_receptor(rigid_pdbqt_filename="bcd_deriv.pdbqt", flex_pdbqt_filename="bcd_deriv_flex.pdbqt")
v.compute_vina_maps(center=[0.,0.,0.], box_size=[25, 25, 25])
v.set_ligand_from_file("hex.pdbqt")
v.dock(exhaustiveness=2, n_poses=1)
scores = v.score()
print(scores[0])
print(v.info())
v.write_poses('output_pose.pdbqt', n_poses=1, overwrite=True)

If I choose flex_pdbqt_filename=None I can dock without issue. (Note that I've also changed np.int to np.int64 as a fix for #196)

Would you be willing to brainstorm ways to solve this problem with me? I'd be nice if there were a way to easily make flexible any sidechain in any structure, just like there's currently a way to easily make any residue flexible in any protein using example/autodock_scripts/prepare_flexreceptor.py.

Best, Jake Anderson

[diogomart]

Hi @jaketanderson

I added an example script here: https://github.com/forlilab/Meeko/tree/reactive/example/cyclodextrin_flexible

It is in branch "reactive" that will be merged soon and released as v0.5. This will change the API and I wanted to write the script for the new API. To install it you'll need to clone the repository, switch to the "reactive" branch, and install from the local repository, e.g. with pip install -e .

The run.py script loads .mol files instead of PDB because PDBs often don't have bond orders, but for the PDB files you posted RDKit seemed to be reading the connect records and keeping correct bond orders.

RDkit has functions to add bonds between molecules, so it is possible to merge the flexible sidechain back with the cyclodextrin core, if needed.

Hope this helps

[jaketanderson]

Thank you so much for putting this together @diogomart! In the next week or so I will test it out and adopt it into my workflow 🙂