Quantifying ternary complex stability through path-based binding free energy calculations

This repository contains the code and inputs required to reproduce the CRBN- and DCAF16-based ternary complex binding free energies reported in:

''Quantifying Ternary Complex Stability Through Binding Free Energy Calculations in the Geometric Route''

Benedict Tan, Alessio Ciulli, Julien Michel

ChemRxiv preprint

All US sample files used to compute the results presented in the paper are also included. The figures presented in the paper can be generated by using the analysis.ipnyb notebook. A full set of functions for analysis (WHAM implementation, equilibration detection using RED, etc.) can be found in under analysis_scripts.

The dependencies for running all python scripts are as follows:

• openmm version 7.7
• numpy
• pandas
• pyyaml
• haddock3
• red-molsim version 0.1

It is also necessary to install the WHAM implementation of Grossfield Labs, and add the wham executable to your PATH.

Note that this repository is quite large and may require a few minutes to clone.

The following files are included for each system:

Unrestrained_MD

• Input structures used for the unrestrained MD simulations
• Distribution of Boresch angles and dihedrals (stored as pickled dictionaries, keys are 'Frames', 'Time (ns)' and 'DOF values')
• runMD.py is the script used to run the unrestrained MD simulation

MMGBSA Calculations

• The input and output files for the MM/GBSA calculations run using Ambertools are stored in the mm_gbsa directory.

All calculations are submitted from the command line as follows (n_rec is the final index of the receptor, for example n_rec=173 for the DCAF16 systems). For the DCAF16-based systems, the strip mask was modified to remove the DDB1 scaffold from the calculation.

ante-MMPBSA.py -p complex.prmtop -c com.prmtop -r rec.prmtop -l ligand.prmtop -s :WAT,CL,NA -n :1-n_rec --radii mbondi2
MMPBSA.py -O -i mmgbsa.in -o results -sp complex.prmtop -cp com.prmtop -rp rec.prmtop -lp ligand.prmtop -y prod.nc

HADDOCK scoring

For each frame extracted from the unrestrained MD simulation, a docking score was generated using the haddock3-score module in HADDOCK3.

US (Umbrella Sampling)

• US_config.yaml contains the specifications of the orientational and conformational restraint Collective Variables (Boresch non-bonded anchor points + equilibrium values, ligand and receptor indices).
• The run_window.py scripts can be used to submit individual US window simulations.
• The separation/SMD.py script can be used to generate the input configurations for the separation US windows.
• The results directories contain the invidual CV samples that are combined using WHAM to obtain the PMFs for each stage of the thermodynamic cycle. The WHAM input metafiles and output files containing the PMF are also included. The script analysis_scripts/us_analysis.py contains functions to rerun WHAM.
• Each system contains the overall binding free energy in the file deltaG.txt, with the breakdown of the individual contributions being listed in deltaG_stages.csv.

Individual RMSD US simulations

The individual stages involved in the bound and bulk state RMSD US simulations are shown below. The integer 'Restrain order' specifies the order in which restraints are turned on, for example in CRBN_len_CK1a, the US simulations for 'CRBN_len_only' correspond to turning on RMSD restraints in lenalidomide-bound CRBN in complex with unrestrained CK1a. In the 'CK1awithCRBN_len' stage, RMSD restraints in CK1a are turned on, with CRBN-lenalidomide already restrained.

BD1_DCAF16

State	Restrain order	RMSD stage name
Bound	0	DCAF16_only
Bound	1	BD1withDCAF16
Bulk	0	DCAF16
Bulk	0	BD1

BD2_DCAF16

State	Restrain order	RMSD stage name
Bound	0	DCAF16_only
Bound	1	BD2withDCAF16
Bulk	0	DCAF16
Bulk	0	BD2

BRD4_DCAF16 and BRD3_DCAF16

State	Restrain order	RMSD stage name
Bound	0	DCAF16_only
Bound	1	BD1withDCAF16
Bound	2	BD2withDCAF16andBD1
Bulk	0	DCAF16
Bulk	0	BD1_only_bulk
Bulk	1	BD2withBD1_bulk

All ternary complexes

State	Restrain order	RMSD stage name
Bound	0	DCAF16_only
Bound	1	BD1withDCAF16
Bound	2	BD2_IBG1withDCAF16andBD1
Bulk	0	DCAF16
Bulk	0	BD1_only_bulk
Bulk	1	BD2_IBG1withBD1_bulk

(Equivalent IBG3 naming for BRD4_IBG3_DCAF16 system)