%%html
<script src="https://bits.csb.pitt.edu/preamble.js"></script>
<style>:root {--jp-cell-prompt-width: 32px}</style>
Structure-Based Drug Discovery¶
Given a protein structure, how do we find a molecule that binds?
Where might it bind?
Identify potential binding pockets.Which molecules should we consider?
Select compounds with suitable shapes and interaction features.How do we rank the candidates?
Score binding poses and prioritize compounds for testing.
Pockets¶
A pocket is a cavity in the protein where a small molecule might bind. Static pockets are apparent in a reference apo structure. Cryptic pockets are not apparent in that structure, but can open transiently during normal protein dynamics and may be stabilized by a bound ligand.
A pocket is considered druggable if it has characteristics similar to known binding sites of small molecules.
Pocket Identification¶
There are a number of tools that use geometric and chemical descriptors to identify pockets given a static structure.
fpocket - works on MD trajectories as well, includes druggability score- CavityPlus - has web interface with additional features
- pyKVFinder - Python module
- DoGSiteScorer - Another web service
fpocket Algorithm¶
- Perform a Voronoi tessellation of heavy atoms
- Identify alpha spheres centered at Voronoi vertices that connect four atoms without containing any.
- Filter and cluster alpha spheres to identify unique cavities
- Small radius in tightly packed regions
- Large radius at flat/convex surface
- Characterize properties of pocket and filter
CAVITY(Plus) Algorithm¶
- Mark 0.5Å grid as occupied or vacant
- Roll a large (10Å radius) sphere around the protein "erasing" vacant points
- Cluster and characterize remaining regions.
SBDD Funnel¶
Pharmacophore¶
IUPAC: The ensemble of steric and electronic features that is necessary to ensure the optimal supra-molecular interactions with a specific biological target structure and to trigger (or to block) its biological response.
- aromatic ring
- hydrophobic area
- positive ionizable
- negative ionizable
- hydrogen bond donor
- hydrogen bond acceptor
Ligand Only: PharmaGist¶
- Decompose ligands into features
- Pairwise align to selected rigit "pivot" molecule
- Combined pairwise alignments to multiple alignment, identify most important features, and select new pivot. Iterate.
- Cluster results to identify most common pharmacophore.
PharmaGist Results¶
Evaluating Screening Results¶
| Actually active | Actually inactive | |
|---|---|---|
| Selected | True positive (TP) | False positive (FP) |
| Not selected | False negative (FN) | True negative (TN) |
Precision: What fraction of selected compounds are active? $$\mathrm{Precision} = \frac{TP}{TP + FP}$$
Recall (True Positive Rate): What fraction of all active compounds did we select? $$\mathrm{Recall} = \frac{TP}{TP + FN}$$
False Positive Rate (FPR): What fraction of all inactives did we select? $$\frac{FP}{FP+TN}$$
F1: Harmonic mean of precision and recall; high only when both are high. $$F_1 = \frac{2\,\mathrm{Precision}\,\mathrm{Recall}} {\mathrm{Precision}+\mathrm{Recall}} = \frac{2TP}{2TP+FP+FN}$$
Receptor Only: PharmMaker (Mixed solvent MD)¶
Probe Simulations¶
Explore potential interactions by simulating small molecular probes.
PyRod, WaterMap: Water is the probe.
Receptor Only: FTMap¶
Exhaustive exploration of probe-protein interactions using FFT followed by clustering.
Recall: 2D Fourier Transform¶
$$\Large X_{u,v} = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{- j 2\pi \left( \frac{ux}{M}+\frac{vy}{N} \right)}$$
Decomposing our function into periodic functions in both dimensions.
from PIL import Image
img = Image.open('imgs/ollycouch.tiff')
img
import scipy, numpy as np, matplotlib.pyplot as plt
import scipy.fft
dogFFT = scipy.fft.fft2(np.array(img))
vals = np.log(np.absolute(dogFFT))
R = vals.shape[0]//2; C = vals.shape[1]//2;
vals = np.roll(vals, shift=R, axis=0); vals = np.roll(vals, shift=C, axis=1)
plt.matshow(vals); plt.xticks([0,C,2*C],labels=[-C,0,C]); plt.yticks([0,R,2*R],labels=[-R,0,R]); plt.colorbar();
Convolutions¶
The integral of the product of two functions as one is shifted.
$$ (f * g)(t) := \int_{-\infty}^\infty f(\tau) g(t - \tau) \, d\tau$$
Convolution Theorem¶
Given two identically sized images $M$ and $S$ (uniformly sampled functions in N-D space):
$$\Large \mathcal{F}(M * S) = \mathcal{F}(M) \cdot \mathcal{F}(S)$$
That is, the Fourier transform ($\mathcal{F}$) of the convolution is equal to the element-wise product of the transforms of the two inputs. Alternatively:
$$\Large M*S = \mathcal{F}^{-1}(\mathcal{F}(M) \cdot \mathcal{F}(S))$$
f = np.zeros((800,800))
f[:10,:10] = 1/100
f_fft = scipy.fft.fft2(f)
vals = np.absolute(f_fft)
R = vals.shape[0]//2; C = vals.shape[1]//2;
vals = np.roll(vals, shift=R, axis=0); vals = np.roll(vals, shift=C, axis=1)
plt.matshow(vals); plt.xticks([0,C,2*C],labels=[-C,0,C]); plt.yticks([0,R,2*R],labels=[-R,0,R]); plt.colorbar();
plt.matshow(np.real(scipy.fft.ifft2(f_fft*dogFFT)),cmap='gray',vmin=0,vmax=255);
Why does this matter?¶
The Fast Fourier Transform algorithm lets us compute the Fourier Transform in $O(N \log(N))$ instead of $O(N^2)$.
$$\Large X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-\frac {j 2\pi}{N}kn}$$
This means we can convolve two 3D grids with complexity $O(N^3\log(N))$ instead of $O(N^6).$
FFT Rigid Body Docking¶
We have two proteins, a receptor (R) and a ligand (L). We want to evaluate every possible translation and rotation of L with respect to R to identify the lowest energy pose.
If we can represent our energy function as a convolution, we can exhaustively compute all (discretized) translations in $O(N^3\log(N))$, which is a tractable calculation.
Simplified Coulomb Convolution:¶
$$\begin{eqnarray} R_p(l,m,n) &=& \sum_i^{N_r} \frac{q_i}{\phi_R(i,l,m,n)} \\ L_p(l,m,n) &=& \begin{cases} q_j &\text{if (l,m,n) overlaps atom with partial charge } q_j \\ 0 & \text{otherwise} \end{cases} \end{eqnarray} $$
Still need to iterate over all rotations.¶
Generalized Fourier methods can also accelerate searches over rotational coordinates
Ligand-Receptor¶
Pharmacophore Search¶
Advantages¶
- Uses expert human insight to define query
- Fast (millions of molecules in seconds)
- Results are already "docked" - aligned to query
Disadvantages¶
- Relies on expert human insight to define query
- Especially difficult when no bound ligand
- Exploration of alternative binding modes requires multiple queries
- Results depend on the conformers and chemical structures in the screened library
- Need a scoring function to productively rank results
Pharmit: Interactive Exploration of Chemical Space¶
Molecular Docking¶
Predict the most likely conformation and pose of a ligand in a protein binding site.
- Sample conformational space
- Score poses
- Ideally score equals affinity or can be used to productively rank compounds
- Score $\ne$ Free Energy
%%html
<iframe width="560" height="315" src="https://3dmol.org/tests/webpages/docking.html" title="docking" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
<script>
$(".jp-InputArea .o:contains(html)").closest('.jp-InputArea').hide();
</script>
Common Approximations in Docking¶
High-throughput docking trades physical detail for speed.
- Receptor flexibility: usually rigid or limited to selected side chains
- Receptor ensembles can represent alternative conformations
- Ligand flexibility: usually samples torsions while fixing bond lengths, bond angles, and ring conformations
- Solvent: usually treated implicitly or through empirical scoring terms
- Selected binding-site waters may be included
- Chemical state: protonation, tautomers, and metal/cofactor treatment depend on input preparation
More detailed treatments are possible, but increase computational cost.
Docking searches for the best scoring pose.
The best score does not guarantee the correct pose or strongest binding.
Protein-Ligand Scoring¶
Force-field based¶
- inter- and intra- molecular forces
- van der Waals, electrostatic, torsional
Empirical¶
- parameterized function is fit to binding energy data
Knowledge based¶
- scoring function based on known structures, not physical principles
Machine Learning¶
- fancy knowledge-based
Empirical: AutoDock Vina Scoring¶
There is no electrostatic term. Partial charges are not used. Electrostatic interactions are (partially) accounted for with hydrogen bond term.
Metals are modeled as hydrogen donors.
Terms were selected and parameterized for pose prediction performance (both speed and quality).
Final scoring function was then linearly reweighted to fit score to free energies (kcal/mol).
Machine Learning: CNN Scoring (GNINA)¶
Atoms are represented as Gaussian densities on a 24Å grid. There is a separate channel for each atom type.
GNINA Docking Algorithm¶
Defining the Binding Site¶
All poses are sampled within a box defined by the user.
Can be specified manually (--center_x, --size_x, etc.) but typically much easier to provide an autobox ligand.
A box is created that exactly encloses the atom coordinates of the provided ligand and then is expanded by autobox_add (default 4Å) in every dimension. If needed to provide enough room for the ligand to freely rotate the box is then further extended (autobox_extend).
Can provide any molecule for autobox_ligand (e.g. binding pocket residues, fpocket alpha spheres).
Sampling¶
Degrees of freedom
- Translation and rotation: 6 rigid-body degrees of freedom
- Rotatable-bond torsions; bond lengths, bond angles, and ring geometry remain fixed
Initialization
- Randomize position, orientation, and torsions
- Fixed geometry, such as ring pucker, still depends on the input conformer
Monte Carlo with local minimization
- Randomly change a translation, rotation, or torsion
- Refine the pose using truncated BFGS with softened scoring potentials
- Apply the Metropolis acceptance criterion
- Always accept a lower score
- Sometimes accept a higher score to help escape local minima
- Retain the best-scoring poses encountered
Scoring¶
Empirical (e.g. Vina)
- Fast and interpretable
- A collection of weighted terms
- By default used for search and refinement
CNN
- Slower (especially without GPU) but more predictive
- By default used only for final ranking
--cnn_scoring=rescore
CNN Scoring Performance¶
Virtual Screening GNINA Performance¶
 |  |
Evaluating Virtual Screening: ROC and EF1%¶
ROC (Receiver Operating Characteristic) curve: plot TPR against FPR as the selection threshold varies.
- ROC AUC summarizes ranking across thresholds: 0.5 for random ranking, 1 for perfect ranking.
EF1% (Enrichment Factor at 1%): how much more concentrated are actives among the top-ranked 1% than in the full library?
$$EF_{1\%} = \frac{\text{active fraction among the top-ranked 1%}} {\text{active fraction in the full library}}$$
ROC AUC measures overall ranking; EF1% emphasizes early retrieval. Random ranking has expected EF1% = 1.
DiffDock Results¶
Example: GDAP1¶
Ganglioside-induced differentiation-associated protein 1 (GDAP1) may play a role in a signal transduction pathway during neuronal development. Mutations in this gene have been associated with various forms of Charcot-Marie-Tooth Disease and neuropathy.
Constrained Cofolding (Boltz)¶
sequences:
- protein:
id: G
sequence: MAERQEEQRGSPPLRAEGKADAEVKLILYHWTHSFSSQKVRLVIAEKALKCEEHDVSLPLSEHNEPWFMRLNSTGEVPVLIHGENIICEATQIIDYLEQTFLDERTPRLMPDKESMYYPRVQHYRELLDSLPMDAYTHGCILHPELTVDSMIPAYATTRIRSQIGNTESELKKLAEENPDLQEAYIAKQKRLKSKLLDHDNVKYLKKILDELEKVLDQVETELQRRNEETPEEGQQPWLCGESFTLADVSLAVTLHRLKFLGFARRNWGNGKRPNLETYYERVLKRKTFNKVLGHVNNILISAVLPTAFRVAKKRAPKVLGTTLVVGLLAGVGYFAFMLFRKRLGSMILAFRPRPNYF # only for protein, dna, rna
- ligand:
id: T # multiple ids in case of multiple identical entities
smiles: 'C(CCCCCCCC(=O)O)CCCCCCC(=O)O'
constraints:
- pocket:
binder: T
contacts: [[G, 36], [G, 134], [G, 257]]
max_distance: 10
force: true # apply a potential to enforce the constraint
templates:
- cif: 6uih.cif # if only a path is provided, will find the best matchings
properties:
- affinity:
binder: T
Boltz Output¶
| confidence_*_.json | affinity_*.json |
|---|---|
{
"confidence_score": 0.865314781665802,
"ptm": 0.8453083634376526,
"iptm": 0.9104629158973694,
"ligand_iptm": 0.9104629158973694,
"protein_iptm": 0.0,
"complex_plddt": 0.8540277481079102,
"complex_iplddt": 0.7315892577171326,
"complex_pde": 0.6098760962486267,
"complex_ipde": 1.484511137008667,
"chains_ptm": {
"0": 0.8421840667724609,
"1": 0.7693953514099121
},
"pair_chains_iptm": {
"0": {
"0": 0.8421840667724609,
"1": 0.39824172854423523
},
"1": {
"0": 0.9104629158973694,
"1": 0.7693953514099121
}
}
}
|
{
"affinity_pred_value": 0.633259654045105,
"affinity_probability_binary": 0.7235137224197388,
"affinity_pred_value1": 0.9157558679580688,
"affinity_probability_binary1": 0.7583574652671814,
"affinity_pred_value2": 0.3507634401321411,
"affinity_probability_binary2": 0.6886700391769409
}
|
Key Points¶
Pockets: identify candidate binding sites; pocket detection alone does not establish druggability.
Pharmacophores: search for ligands matching a 3D arrangement of interaction features.
Docking = sampling + scoring: search ligand positions, orientations, and conformations.
- FFT accelerates evaluation of rigid-body translations.
Scoring: force-field, empirical, statistical, and machine-learning approaches.
- A favorable score or high confidence does not establish binding.
Evaluation: pose accuracy, affinity prediction, and virtual screening are different tasks.
- ROC AUC measures overall ranking; EF1% measures early enrichment.
Deep learning: improves scoring and enables learned docking and cofolding.
- Assess generalization to new targets and chemotypes, accounting for training overlap.